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Abstract 

The Booster Neutrino Experiment (MiniBooNE) searches for — > v e oscillations using the 
(9(1 GeV) neutrino beam produced by the Booster synchrotron at the Fermi National Accelerator 
Laboratory (FNAL). The Booster delivers protons with 8 GeV kinetic energy (8.89 GeV/c momen- 
tum) to a beryllium target, producing neutrinos from the decay of secondary particles in the beam 
line. We describe the Monte Carlo simulation methods used to estimate the flux of neutrinos from 
the beamline incident on the MiniBooNE detector for both polarities of the focusing horn. The 
simulation uses the Geant4 framework for propagating particles, accounting for electromagnetic 
processes and hadronic interactions in the beamline materials, as well as the decay of particles. 
The absolute double differential cross sections of pion and kaon production in the simulation have 
been tuned to match external measurements, as have the hadronic cross sections for nucleons and 
pions. The statistical precision of the flux predictions is enhanced through reweighting and re- 
sampling techniques. Systematic errors in the flux estimation have been determined by varying 
parameters within their uncertainties, accounting for correlations where appropriate. 
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I. INTRODUCTION 



The Booster Neutrino Experiment (MiniBooNE) at Fermilab searches for the oscillation 
of muon neutrinos (i/^) to electron neutrinos (f e ) indicated by the LSND experiment [1J[2J. 
The neutrino beam is produced by the Booster Neutrino Beamline (BNB), where protons 
with 8 GeV kinetic energy (8.89 GeV/c momentum) are extracted from the Fermilab Booster 
synchrotron and directed towards a beryllium target. Secondary mesons produced by the 
interaction of the protons in the target decay to produce a neutrino beam with an average 
energy of ~ 800 MeV. Neutrino interactions are observed in a 6.1-meter-radius, spherical 
detector situated 541 meters from the center of the target. The detector is composed of 
800 metric tons of mineral oil that serves as both the target for neutrino interactions and 
the medium in which charged particles produced in neutrino interactions radiate Cherenkov 
and scintillation photons. The photons are detected on an array of 1520 photomultipliers, 
and the resulting spatial and temporal patterns of light are used to identify and reconstruct 
the interactions. Understanding both the spectrum and composition of the neutrino beam 
is critical to the neutrino oscillation analysis, which searches for an excess of u e events over 
a background of both non-oscillation sources of u e in the beamline and misidentified 
interactions. 

The neutrino oscillation analysis at MiniBooNE utilizes observed data to constrain the 
uncertainties in the expected event rates of certain key processes. These constraints typically 
reduce the uncertainties that would result from a direct estimation using solely the predicted 
neutrino flux and cross sections. Within the — > u e oscillation analysis, the observed rate 
of charged current quasi-elastic events and neutral current tt° events are used directly in 
the estimation of the number and spectrum of background and expected neutrino oscillation 
events. The Monte Carlo-based flux prediction described here is one input to this process. 
In this article, we focus on the flux prediction itself, which is based on external data, without 
regard to the observed neutrino rates at MiniBooNE. Predictions for both polarities of the 
focusing horn are presented. Detailed comparisons of the observed and predicted neutrino 
event rates, as well as descriptions of the use of the predicted neutrino fluxes in various 
analyses, are described in publications relating to the analysis of the neutrino data itself, 
e.g., References [21 El H [5]. 
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II. THE BOOSTER NEUTRINO BEAMLINE 



The Booster Neutrino Beamline (BNB) produces neutrinos using 8.89GeV/c momentum 
protons from the Booster synchrotron that are incident on a beryllium target. The layout of 
the BNB is shown in Figure [TJ The target is embedded within a pulsed electromagnet (the 
"horn") that produces a toroidal magnetic field to focus positive secondary particles and 
defocus negative secondary particles emerging from proton-beryllium interactions. These 
secondary particles enter a 50-meter-long decay region, resulting in a neutrino-enhanced 
beam. The polarity of the horn can be reversed to focus negative secondary particles and 
produce an antineutrino-enhanced beam. The axis of the beam, defined by the center of the 
decay pipe, is displaced vertically from the center of the MiniBooNE detector by 1.9 meters. 

The particle production is dominated by pions, though there is significant kaon produc- 
tion as well. Neutrinos also result from the decay of muons whose primary source is the 
decay of pions produced in the target. This results in a significant flux of z/ e /z7 e in neu- 
trino/antineutrino mode, while the corresponding flux of v^jv^ is small compared to the 
v iil^ii which result directly from the decay of the pions. A beam stop at the end of the 
decay region absorbs particles apart from the neutrinos. The predicted composition of the 
neutrino beam is described in Section fVIj A detailed description of the BNB can be found in 
the Technical Design Report for the BNB [6]. This section describes the beamline geometry 
and components relevant for the neutrino flux prediction. 

A. FNAL Booster and Proton Extraction 

The FNAL Booster is a 474-meter-circumference synchrotron operating at 15Hz. Protons 
from the Fermilab LINAC are injected at 400 MeV and accelerated to 8 GeV kinetic energy 
(8.89GeV/c momentum). The Booster has a harmonic number of 84, of which 81 buckets 
are filled. The beam is extracted into the BNB using a fast-rising kicker that extracts all 
of the particles in a single turn. The resulting structure is a series of 81 bunches of protons 
each ~ 2 ns wide and 19 ns apart. 

Upon leaving the Booster, the proton beam is transported through a lattice of focusing 
and defocusing quadrupole (FODO) and dipole magnets. A switch magnet steers the beam 
to the Main Injector or to the BNB. The BNB is also a FODO that terminates with a 
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FIG. 1: Overall layout of the BNB. The primary proton beam, extracted from the Booster, enters 
the target hall from the left. Upon exiting the target hall, particles encounter a 50-meter-long 
decay region, terminating in the beam stop on the right. 
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FIG. 2: The MiniBooNE target assembly. The top shows an exploded view of the components, 
while the bottom shows the assembled configuration. The proton beam enters from the left in both 
figures, striking the finned beryllium slugs. Dimensions are in inches. 



triplet that focuses the beam on the target. The design and measured optics of BNB are in 
agreement [H] . 

The maximum allowable average repetition rate for delivery of protons to the BNB is 5 
Hz (with a maximum of 11 pulses in a row at 15 Hz) and 5 x 10 12 protons-per-pulse. The 
5 Hz limit is set by the design of the horn (described below) and its power supply. As of 
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FIG. 3: Left: Neutrino event times relative to the nearest RF bucket (measured by the RWM) 
corrected for expected time-of- flight. Right: An oscilloscope trace showing the coincidence of the 
beam delivery with the horn pulse. The top trace (labeled "2" on the left) is a discriminated signal 
from the resistive wall monitor (RWM), indicating the arrival of the beam pulse. The bottom trace 
(labeled "1" on the left) is the horn pulse. The horizontal divisions are 20 fjs each. 

January 2008, over 10 21 protons have been delivered to the BNB, with a typical up time 
of greater than 90% during normal operations. The neutrino oscillation results in neutrino 
mode were published using 5.6 x 10 20 protons-on-target delivered prior to 2006, when the 
polarity of the horn was reversed to collect antineutrino mode data [9]. 

B. Target 

The target consists of seven identical cylindrical slugs of beryllium arranged to produce 
a cylinder 71.1 cm long and 0.51 cm in radius. The target is contained within a beryllium 
sleeve 0.9 cm thick with an inner radius of 1.37 cm. Each target slug is supported within 
the sleeve by three "fins" (also beryllium) which extend radially out from the target to 
the sleeve. The volume of air within the sleeve is circulated to provide cooling for the 
target when the beamline is in operation. The target and associated assembly are shown in 
Figure [2j where the top figure shows an "exploded" view of the various components (with 
the downstream end of the target on the right), and the bottom shows the components in 
assembled form. The choice of beryllium as the target material was motivated by residual 
radioactivity issues in the event that the target assembly needed to be replaced, as well as 
energy loss considerations that allow the air-cooling system to be sufficient. 



7 



Upstream of the target, the primary proton beam is monitored using four systems: two 
toroids measuring its intensity (protons-per-pulse), beam position monitors (BPM) and a 
multi-wire chamber determining the beam width and position, and a resistive wall monitor 
(RWM) measuring both the time and intensity of the beam spills. The vacuum of the 
beam pipe extends to about 5 feet upstream of the target, minimizing upstream proton 
interactions. 

The toroids are continuously calibrated at 5 Hz with their absolute calibrations verified 
twice a year. The calibrations have shown minimal deviation (< 0.5%). The proton flux 
measured in the two toroids agree to within 2%, compatible with the expected systematic 
uncertainties. The BPMs are split-plate devices that measure the difference of charge in- 
duced on two plates. By measuring the change in beam position at several locations without 
intervening optics, the BPMs are found to be accurate to 0.1 mm (standard deviation). The 
multi-wire is a wire chamber with 48 horizontal and 48 vertical wires and 0.5 mm pitch. 
The profile of the beam is measured using the secondary emission induced by the beam on 
the wires. 

The RWM is located upstream of the target to monitor the time and intensity of the 
proton pulses prior to striking the target. While the data from the RWM did not directly 
enter the — > v e analysis, it allowed many useful cross checks, such as those shown in 
Figure [3] The left figure shows a comparison of the production times of neutrinos observed 
in the MiniBooNE detector estimated based on the vertex and time reconstructed by the 
detector and subtracting the time-of-flight. This time is then compared to the nearest bucket 
as measured by the RWM. The distribution indicates that neutrino events can be matched 
not only to pulses from the Booster, but to a specific bucket within the pulse. The tails 
of the distribution result from the resolution of the vertex reconstruction of the neutrino 
event in the detector, which is needed to determine the time of the event and correct for the 
time-of-flight. The right plot demonstrates the synchronization of the horn pulse (described 



in Section II C) with the delivery of the beam as measured by the RWM. 



C. Horn Electromagnet 



The horn, shown in Figure |4j is a pulsed toroidal electromagnet composed of an aluminum 
alloy (6061 T6). The current in the horn is a 143 //s-long pulse with a nominal peak of 170 
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FIG. 4: The MiniBooNE pulsed horn system. The outer conductor (gray) is transparent to show the 

inner conductor components running along the center (dark green and blue). The target assembly 

is inserted into the inner conductor from the left side. In neutrino-focusing mode, the (positive) 

current flows from left-to-right along the inner conductor, returning along the outer conductor. 

The plumbing associated with the water cooling system is also shown. 
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FIG. 5: Measurements of the azimuthal magnetic field within the horn. The points show the 
measured magnetic field, while the line shows the expected 1/R dependence. The black lines 
indicate the minimum and maximum radii of the inner conductor. 



kA coinciding with the arrival of the proton beam at the target. The actual operating values 
are typically 174 kA in both neutrino and antineutrino mode with ~ 1 kA variations. In 
neutrino mode, the flow of current (in the positive sense) runs along the inner conductor 
(containing the target), which folds outwards at a length of 185 cm to return via the the 
outer conducting cylinder of the horn at 30 cm radius. Within the horn cavity, defined by the 
volume between the outer and inner conducting cylinders, the pulse creates a magnetic field 
that falls as 1/R, where R is the distance from the cylindrical symmetry axis of the horn. 
The largest field values of 1.5 Tesla are obtained where the inner conductor is narrowest 
(2.2 cm radius). The effects of time- varying fields within the cavity of the horn are found 
to be negligible. The expected field properties of the horn have been verified by measuring 
the current induced in a wire coil inserted into the portals of the horn. Figure [5] shows the 
measured R dependence of the azimuthal magnetic field compared with the expected 1/R 
dependence. The "skin effect", in which the time- varying currents traveling on the surface 
of the conductor penetrate into the conductor, results in electromagnetic fields within the 
conductor itself. 

During operation, the horn is cooled by a closed water system which sprays water onto 
the inner conductor via portholes in the outer cylinder. The target assembly is rigidly fixed 
to the upstream face of the horn, although the target is electrically isolated from its current 
path. At the time of writing, two horns have been in operation in the BNB. The first 
operated for 96 million pulses before failing, while the second is still in operation as of this 
writing after over 130 million pulses [TU] . 

D. Collimator 

A concrete collimator is located downstream of the target/horn assembly and serves as the 
entrance into the decay region. The collimator absorbs particles that would not otherwise 
contribute to the neutrino flux and is 214 cm long, with an upstream aperture of 30 cm radius 
that grows to 35.5 cm on the downstream end. By absorbing these particles, the collimator 
reduces radiation elsewhere in the beamline.The upstream end of the collimator is located 
259 cm from the upstream face of the target. Simulations indicate that the collimator does 
not limit the neutrino flux at the MiniBooNE detector. 
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E. Decay Region and Absorber 

The air-filled cylindrical decay region following the collimator is 3 feet in radius and 
extends for 45 meters, terminated by the beam stop 50 meters from the upstream face of the 
target. Survey data indicate that the constructed decay region is 49.87 meters, including 
the collimator region. The wall of the decay region is a corrugated steel pipe surrounded 
by packed dolomite gravel (CaMg(C0 2 ) 3 , p = 2.24g/ cm 3 ). The beam stop itself is made of 
steel and concrete, within which is an array of gas proportional counters that detects muons 
penetrating the beam stop. 

To allow potential systematic studies, a set of ten steel absorbing plates are positioned 
above the decay pipe at 25 meters. When lowered into the decay region, the steel absorbers 
reduce the effective decay path from 50 to 25 meters. This has the effect of reducing the 
overall flux, but preferentially reducing the decay of the longer-lived muons, a major source of 
non-oscillation v e background. The 25 meter absorber was not deployed during the neutrino 
running for the — > v e oscillation analysis. 

F. Little Muon Counter 

The Little Muon Counter (LMC) is an off-axis spectrometer that measures the rate and 
spectrum of muons produced at a 7° angle to the beam axis in the decay pipe pointing back to 
the alcove for the 25 meter absorber. The detector consists of a 40 foot drift pipe extending 
from the decay region at 7° leading to an enclosure. The kinematics of the two-body pion and 
kaon decay are such that kaons produce a momentum distribution peaked at 1.8 GeV/c at this 
angle, whereas pions produce muons at lower momentum. Muons sent down the drift pipe 
to the enclosure encounter an iron collimator with a tungsten core that further restricts the 
angular acceptance of the counter and reduces backgrounds. Following the collimator, the 
muon momentum is determined by a spectrometer consisting of a dipole magnet and planes 
of scintillating fiber trackers. Finally, a range stack consisting of alternating scintillator and 
tungsten layers allow high energy muons to be distinguished from pions and other particles 
based on the number of tungsten planes penetrated by the particle. Further details on the 
LMC can be found in Reference |llj . 
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III. GEANT4-BASED MONTE CARLO SIMULATION 



The properties of the MiniBooNE neutrino flux are determined using a Geant4-based 
Monte Carlo simulation [121 . The simulation can be divided roughly into five steps: 



The definition of the beamline geometry, specified by the shape, location, and com- 
position of the components of the BNB, through which the primary protons and all 



other particles propagate (Section III A). 



The generation of the primary protons according to the expected beam optics proper- 



ties upstream of the target (Section III C ) . 



The simulation of particles produced in the primary p-Be interactions, including the 
elastic and quasi-elastic scattering of protons in the target. Custom tables for the 
production of protons, neutrons, 7r , K and K° in these interactions have been 



developed to accommodate production models based on external data (Section HIE) 



The propagation of the particles using the Geant4 framework, accounting for energy 
loss and electromagnetic and hadronic processes that alter the kinematics of the parti- 



cles as described in Section HI D Hadronic interactions and decay processes may also 
annihilate the particle in the tracking process and create new particles to be tracked. 
Within the horn, the effect of the expected magnetic field on the trajectory of the 



particles is accounted for (Section III B ) 



The identification of decay processes that result in neutrinos. The simulation of the 



decays is handled by a custom decay model, described in Section |III F[ outside of the 
Geant4 framework. The decay model reflects the latest branching fraction measure- 
ments and simulates polarization effects and kinematic distributions resulting from 
decay form factors. A number of techniques to enhance the statistical precision of the 



flux prediction are employed (Section IIIG). 



A. Geant4 Description of the BNB Geometry: 

The Geant4 Monte Carlo geometry consists of the last 50 meters of the Booster beam- 
line, the target hall, and the 50 meter meson decay volume. The geometry description is 
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FIG. 6: The MiniBooNE horn as simulated in the Geant4 beam Monte Carlo. 



defined to match the actual constructed beamline as closely as possible; differences from 
the specifications are noted here. Since the generation of primary protons in the simulation 



starts immediately upstream of the target (see Section III C ) , the geometry description of 
the beamline leading to the target is simplified. Each section is simulated with concrete 
walls surrounded by a uniform bed of dolomite. The entire structure is filled with air at 
standard temperature and pressure. 

The simulated target hall contains the target, horn, and secondary beam collimator. In 
addition to the concrete walls, the target hall is lined with 1.28 meters of steel shielding. 
The seven slugs of the target and the target sleeve, together with the fins which support 
them within the sleeve, compose the simulated target geometry. The horn is constructed 
using an aluminum Geant polycone that specifies the inner and outer radius at 14 different 
points along the direction of the beam. A polycone of air is placed inside of the aluminum 
polycone to set the thickness of the inner and outer conductors. A graphical representation 
of the Monte Carlo horn is shown in Figure [6j 

The meson decay region is simulated as two 20-meter-long decay pipes separated by the 
25 meter absorber enclosure, followed by the beam stop. The decay pipes are made of 
concrete with an inner diameter of 6 feet and an outer diameter of 10 feet. In contrast, the 
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actual decay pipe is corrugated steel and surrounded by dolomite; this simplification is not 
expected to affect the flux prediction. 



B. Horn Magnetic Field 

The 1/R magnetic field expected from the current path within the horn is simulated in 
the volume corresponding to the cavity of the horn. The strength of the field corresponds to 
a 174 kA current running along the inner conductor (reversed when simulating anti- neutrino 
mode). In addition, the permeation of the magnetic field into the inner conducting cylinder 



of the horn from the skin effect (described in Section VII ) is included in the simulation. The 
predicted trajectories for charged particles in the magnetic field in the Geant4 simulation 
have been checked in an external study using the DRKNYS routine from CERNLIB [13J, 
an independently implemented numerical method. 

C. Generation of the Primary Proton Beam 

The primary protons are simulated individually, since no correlated effects between the 
protons in a bunch are expected. The properties of the proton beam, such as the position 
and profile, have been simulated using TRANSPORT [H] and verified by upstream beam 
monitors [8j. The protons are generated 1 cm upstream of the target with the transverse 
(x, y) positions drawn from random Gaussian distributions with mm mean and 1.51 mm 
and 0.75 mm widths, respectively. Likewise, the angular deviations of the proton direction 
from the z direction, 9 X and 9 y , are drawn from Gaussian distributions with mrad mean 
and 0.66 mrad and 0.40 mrad width, respectively. The number of protons that undergo 



inelastic interactions in the target (as opposed to scattering out) is studied in Section VII A 
In particular, while the default configuration describes a diverging beam, the TRANSPORT 
simulations indicate that the protons are expected to be convergent on the target, with a 
"waist" of zero divergence at the center of the target. The simulated beam configuration is 
such that 99.8% of the protons are on a trajectory to intersect the target. The studies in 



Section VII A indicate that reasonable perturbations to the model, including the expected 



focusing configuration, do not affect the predicted neutrino flux by more than 1%. 
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D. Particle Tracking and Propagation 



Particles, whether they are primary protons or particles produced in the simulation, are 
tracked and propagated within the Geant4 framework with full accounting for electromag- 
netic and hadronic interactions and decays. Within each medium, the Coulomb scattering 
and energy loss are calculated in each step of the tracking, and the particle trajectory and 
energy are updated accordingly. The energy loss and deflection angles predicted by the 
framework have been checked in a comparison with the Bethe-Bloch formalism and the 
Highland formula [TU] . The rate of hadronic interactions for protons, neutrons and charged 
pions on beryllium and aluminum are governed by customized cross section tables (see Sec- 
tion IV) that determine the rate of elastic and inelastic scattering within the target and 



horn. The outgoing final state configurations of these interactions are handled by the de- 
fault Geant4 elastic and inelastic scattering algorithms, with the exception of the primary 
p-Be interactions. 

For other particle/nucleus combinations, the default cross section tables in Geant4 are 
used. Extensive checks have been performed to ensure that the rate of hadronic interactions, 
both elastic and inelastic, are consistent with cross sections assigned to these processes. The 
final state configurations of neutrinos produced in decays are handled outside of the Geant4 
framework as described in Section IIII Gl 



E. Primary Proton Interactions 

For the vast majority of primary protons, the first beamline component encountered is 
the target. Since p-Be interactions are the primary source of secondary mesons, a dedicated 
model (described in Section [V]) tuned to external data is used to describe the particle 
production in proton interactions in the target slugs, fins and sleeve. Elastic scattering 
is handled by existing Geant4 models, while nucleon and pion quasi-elastic scattering use a 
dedicated model based on free hadron-nucleus scattering data. Due to the divergence of the 
primary beam and scattering, it is possible for a primary proton to interact outside of the 
target (usually the aluminum of the horn or the concrete of the decay region). For these 
cases, the particle production is handled by the default Geant4 hadronic model. 

For the primary p-Be interactions, secondary particles of seven types (7r ± , K ± , K°, p, 
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n) are generated according to custom production tables describing the double differential 
cross sections for the production of each secondary species as a function of p z and p?, the 
components of momentum along and transverse to the primary proton direction, respectively. 
The total production cross section for a given species, obtained by integrating the double 
differential cross section, determines the average multiplicity of the species in each primary 
p-Be reaction (hadronic interactions excluding elastic and quasi-elastic scattering) when 
divided by the reaction cross section. In each such reaction, the multiplicity for each species 
in drawn from a Poisson distribution based on the average multiplicity, and the kinematics 
drawn from the table of double differential cross sections. The reaction cross sections are 
described in Section [IVJ while the specification of the double differential cross section tables 
is described in Section |V1 

F. Particle decays 

Neutrinos reaching the MiniBooNE detector are produced in the decays of charged pions, 
charged and neutral kaons, and muons. The particle lifetimes, decay modes and associ- 
ated branching ratios, and kinematic distributions of the neutrinos produced in the decays 
assumed in the simulation affect neutrino flux predictions, and are discussed here. The 
neutrino parent lifetimes and branching ratios used in the simulation are given in Table. [TJ, 
for 7r + , K + , K®, and fi + neutrino parents. The corresponding decays of negatively-charged 
particles are also simulated. 

In the two-body decays of charged pseudo-scalar mesons M + — ► l + + z//, where 
M = 7r, K and I = e, //, neutrinos are produced in the meson rest frame with fixed en- 
ergy E v = {m 2 M — mf)/ (2rriM) and isotropic angular distribution. 

For kaon semileptonic decays, K — > ir + I + v\ ("Kj 3 "), neutrinos are produced with 
isotropic angular distribution in the kaon rest frame. For the neutrino energy distributions, 
different parametrizations are used for the K l3 and Kf 3 form factors depending on whether 
electron or muon neutrinos are produced in the decay. In both cases, we assume that only 
vector currents contribute, and that time-reversal invariance holds. 

For Kf 3 and K® 3 decays, the electron neutrino energy distribution in the kaon rest frame 
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TABLE I: Particle lifetimes, and neutrino-producing decay modes and branching ratios considered 
in the simulation. 



is given by [151 : 



— cx / dE e (2E e E u - m k E> n )\f e + (t)f 
dE v J E _ 



(1) 



where I£ e is the electron energy, / + is a form factor depending on the square of the four- 
momentum transfer to the pion, t = (j>k — Ptt) 2 = rn 2 K + m l ~ ^ m KE n , E' n is given by: 



E , _ m \ + m l ~ m l _ E 



2rriK 

and E e> ± are integration limits on the electron energy: 

- — 2m t ^ 

We assume a linear dependence of the form factor /i on t: 



(2) 



(3) 



fUt) = f e JO)(l + XU/ml) 



(4) 



For both K^ 3 and _K"g 3 decays, the coefficient J\fj_ for the linear expansion of the form 
factor used is 2.82 ■ 1(T 2 [16]. 
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For K~s and K® 3 decays, the muon neutrino energy distribution in the kaon rest frame 



is given by 



where: 



dN 
dE„ 



oc 



dE,\f^t)\ 2 [A + Bat) + Cat? 



(5) 



A = m K (2E„E v - m^) + ml {E'JA - E v ) 



(6) 



B = ml(E v - E'J2) 
C ee m^;/4 

The quantities £^ _, /+(£) appearing in Equations [5] and [6] are defined as 

in the K e3 case (see Equations [2j |3j substituting e — ► /i. In the simulation, we take 
f(t) ~f(0) = -0.19 US]. 

Concerning electron (anti-)neutrinos from \i decays, we neglect terms proportional to 
the electron mass and assume the following neutrino energy and angular distribution [T7j : 

dN 12x 2 



[1-x){1tP z cos(6 u )] 



(7) 



dxdQy 4:71 

where cos(8 u ) is the neutrino emission angle with respect to the beam direction z, P z is the 
projection along z of the muon polarization vector in the muon rest frame, and x = 2E u /m^, 
with < x < 1. The muon polarization vector is estimated on an event-by-event basis. For 
7r + — > fi + — > z/ e decays, the muon polarization in the muon rest frame is calculated from 
the known muon polarization in the pion rest frame, and boosting the polarization vector 
into the muon rest frame. The muon polarization for muons proceeding from K ± decays is 
computed in the same way, with the simplifying assumption that all K decays proceed via 
the K^ 2 decay mode. 



G. Statistical Enhancements 



Running the Geant4 beamline simulation and recording the outgoing neutrinos proton- 
by-proton would not provide enough neutrinos at the MiniBooNE detector to allow for a 
precise determination of the flux across the entire phase space of interest. As a result, several 
modifications are made to enhance the beam Monte Carlo simulation statistics. 

A large statistical enhancement is gained by "redecaying" the parent particle of the 
neutrino. For each neutrino produced in the beam Monte Carlo, the particle decay which 
produced the neutrino is performed 1000 times. Each redecay is performed at the same 
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location, but the kinematics of the decay are randomly redrawn each time from the decay 
distributions, resulting in different momenta for the daughter particles in each draw. 

A similar technique is used to boost the statistics of neutrinos from muon decay. Most 
muons produced in the secondary beam do not decay before stopping in the beam stop 
or the walls of the decay region due to their long lifetime. To better estimate the muon 
decay-in-flight component of the neutrino flux, each time a muon is produced in the simula- 
tion, 19 identical copies are created and independently propagated through the simulation. 
To account for the resulting overproduction of neutrinos, the weight for each muon-decay 
neutrino is correspondingly reduced by a factor of 20. 

Another weighting technique is used to determine the high energy neutrino flux. The 
MiniBooNE neutrino energy spectrum peaks between 500 and 600 MeV with a long high 
energy tail extending past 6 GeV. Since fewer neutrinos are produced at these higher ener- 
gies, statistical fluctuations are much larger, increasing the uncertainty in the shape of the 
high energy tail. This problem is made worse by the redecay procedure described above 
since high energy parent particles tend to decay to high-energy, forward-going neutrinos, 
which resulting in a significant fraction of the 1000 redecays producing neutrinos pointed 
at the detector, all with similar energies. To reduce the statistical uncertainty in the pre- 
diction of the high-energy tail, the meson production cross sections for proton-beryllium 
interactions are multiplied by an exponential function in longitudinal meson momentum. 
Each event is de-weighted by its corresponding cross section enhancement to preserve the 
correct neutrino/proton ratio. This provides an artificially large production of neutrinos at 
high energies with small event weights, thus reducing the statistical uncertainty in the high 
energy tail of the predicted neutrino flux. 

The neutrino flux at the MiniBooNE detector is determined by projecting the path of 
the neutrino to the plane containing the center of the detector, 541 meters from the face 
of the target. Neutrinos which are on a path to pass through the detector (within 610.6 
cm of the center of the detector at this plane, accounting for the vertical displacement) are 
recorded in the flux distributions used for Monte Carlo event simulation in the detector. 
For the simulation of neutrino interactions outside the detector in the concrete walls of the 
detector hall or in the dirt beyond, a larger radius of 1400 cm is used to determine the flux 
distributions. 



19 



IV. HADRONIC INTERACTIONS IN THE BEAMLINE 



Hadronic cross sections play an important role in determining the properties of the neutri- 
nos produced in the BNB. Most notably, hadronic cross sections on beryllium and aluminum, 
the materials composing the target and horn, respectively, govern the rate of primary proton 
interactions in the target, as well as the rate of absorption of pions in the target and the 
horn. The breakdown of the proton cross sections between elastic, quasi-elastic and other 
forms of interactions govern the fraction of protons that scatter out of the target before 
interacting. The analogous breakdown of the cross section for pions is particularly impor- 
tant at high momentum, where forward-going pions may intersect a considerable amount of 
material in the target or the horn before entering the decay region. 

The cross sections fall into three categories: elastic (coherent) scattering, inelastic scat- 
tering, and quasi-elastic scattering. In the first, the incident hadron scatters coherently 
from the nucleus as a whole. The rest of the total hadronic cross section is due to inelastic 
processes. A subset of these processes involve hadron scattering with the nucleons within 
the nucleus in a manner analogous to the elastic scattering of hadrons off free nucleons; 
this is referred to as quasi-elastic scattering. The remainder of the inelastic cross section 
includes the particle production processes discussed in Section |Vj The relevant momentum 
range in the flux prediction are at and below the primary proton momentum (8.89GeV/c) 
for nucleons. The corresponding momentum range for the pions produced by these protons 
is - 6 GeV/c 

Wherever possible, measured cross sections have been used in the simulation. In some 
cases, the measured and calculated cross sections are extrapolated to cross sections for other 
particles that are related by isospin. Measurements exist primarily for total hadronic cross 
sections and inelastic cross sections, from which the elastic cross section can be inferred. 
Theoretical guidance is needed primarily for the total hadronic cross sections for pion- 
nucleus scattering and quasi-elastic scattering. Table [TT] summarizes the source of nucleon 
and pion cross sections. Details of the parametrizations used to describe the momentum 



dependence of each cross section are given in Section IV C 
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p-(Be/Al) 


n-(Be/Al) 


^-(Be/Al) 


cttot 


Glauber 


Glauber 


Data (p < 0.6/0.8 GeV/c) 






(checked with data) 


Glauber (p > 0.6/0.8 GeV/c) 


vine 


Data 


(same as p-Be/Al) 


Data 


&QEL 


Shadow 


Shadow 


Data (p < 0.5 GeV/c ) 








Shadow (p > 0.5 GeV/c) 



TABLE II: Origin of hadron-beryllium cross sections used in the Geant4 simulation. "Glauber" 
indicates that the Glauber model calculations described in Section IIV Al are used. These have 
been cross checked by n-Be gtot data. "Data" indicates that existing measurements are directly 
parametrized. "Shadow" refers to the calculation of <jqel using the shadowed multiple scattering 
model described in Section IIV Bl 



A. Total and Elastic Scattering 

The elastic scattering cross sections for protons, neutrons and charged pions have been 
obtained by calculating the total hadronic cross section <Jtot using the Glauber model 



and subtracting the measured inelastic cross sections described in Section IV B assuming 
<Jtot = pine + &ela- Direct measurements of otot for hadron-nucleus interactions in the 
relevant energy range exist only for neutrons and for pions in the A (1232) resonance region. 
Wherever possible, we compare the calculated results with the existing measurements to 
check their validity. No direct measurements of the total elastic cross section {(Tela) exist 
in this momentum range. 

The calculation of otot follows the work described in Reference |24j, where hadron- 
nucleus elastic scattering is modeled as the coherent sum of scattering amplitudes from 
hadron-nucleon scattering. The amplitude for forward elastic scattering is calculated al- 
lowing Otot to be obtained via the optical theorem. At each incident hadron energy, these 
amplitudes are summarized by three parameters, namely the total cross section for hadron- 
nucleon scattering (a n ), the ratio of the real and imaginary parts of the forward scattering 
amplitude (a), and the differential cross section in t = \q 2 \, the square of the 4-momentum 
transfer. The latter is parametrized as an exponential distribution in t. All together, the 
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FIG. 7: Total and elastic hadron-proton cross sections (Top: proton/neutron, Bottom: 7r + /7r~) 
compiled by the Particle Data Group and the COMPAS collaboration[18], with the parametriza- 
tions used in the Glauber model calculations. 

hadron-nucleon scattering amplitude can be expressed as: 

= (i + a)ka n ePt/2 
An 

where k is the wave number of the incident hadron. This form identically satisfies the optical 
theorem. 

The Glauber model for elastic scattering represents the nucleus as a collection of nu- 
cleons distributed in a spherically symmetric state with radial distributions given by the 
independent harmonic oscillator form (for beryllium) or the Woods-Saxon form (for alu- 
minum) (25]. The scattering amplitude for a given configuration of nucleons is obtained 
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Momentum (GeV/c) Momentum (GeV/c) 



FIG. 8: Measured values of a the real-to-imaginary ratio of the forward scattering amplitude for 
p — p and n — p scattering (left) and n + — p and 7r~ — p scattering (right) with parametrizations. 
The parametrizations used in the Glauber model calculation are shown. 




Momentum (GeV/c) Momentum (GeV/c) 



FIG. 9: (3 parameters for p — p (left) and n — p (right) scattering obtained from fits to the data 
with the parametrizations used in the Glauber model calculation. 

by considering the phase shift due to the individual hadron-nucleon scattering amplitudes. 
The total scattering amplitude for the nucleus is calculated by averaging over all nucleon 
configurations weighted by the nucleon density distribution. The total cross section <Jtot 
is obtained by applying the optical theorem to the resulting forward scattering amplitude. 
As mentioned above, oela at a particular incident hadron momentum is calculated via the 
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FIG. 10: Compilation of measured j3 parameters for n + — p elastic scattering (left) and ir~ — p 
elastic scattering (right) versus incident pion momentum with the parametrizations used in the 
Glauber model calculation. The measured values of (3 include the Lasinski compilation |19j as well 
as our own fits to the q 2 distributions measured in References |20[ |2T| [22] . 

relation otot = vela + vine using the values of (Tine described below. While it is in princi- 
ple possible to obtain a el a from the Glauber model by obtaining the elastic cross section as 
a function of q 2 and integrating, the assumptions of the model are most valid in the forward 
direction. An extraction of a el a using the model requires integrating the differential cross 
section outside of this region. 

The hadron-nucleon scattering parameters a n , a and (3 are obtained from the literature. 
In particular, a n and a for p — p, p — n, n + — p and 7r~ — p elastic scattering have been 
compiled by the Particle Data Group (PDG) [18J. The compiled data on a n and a and our 
parametrization of their momentum dependences are shown in Figure [7] and [8] In addition 
to the PDG compilation, a compilation of a measurements by CERN|22J, as well as the 
measurements of Foley et al. |26j, have been included in the parametrization. This latter 
data are at momenta above the region of interest ([7— 10] GeV/c) but are nonetheless useful 
in determining the momentum evolution of a for 7r ± — p scattering. The measurements 
at momenta less than 3.5GeV/c come entirely from the CERN compilation. The PDG 
compilation of a in nucleon-nucleon scattering adequately covers the range of interest for 
the flux prediction. 

Unfortunately, the f3 measurements for hadron-nucleon elastic scattering have not been 
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FIG. 11: Total hadronic cross sections calculated using the Glauber model and the optical theorem 
for beryllium (left) and aluminum (right) targets. The calculated results for neutrons(protons) are 
shown as triangles(circles). The parameterization used in the flux prediction is shown as the solid 
line, while the default Geant4 parametrization is shown as a dashed line. The measured values of 
otot for n-Be/Al from References [32j |3H [35j [3B] are shown as squares. 



compiled by the PDG. We have taken data from a number of experiments (for p — p[2T] [271 
23 ES], for n - p [301 ED], and for tt ± - p [H EH E2J) and used the reported t distributions 
to extract (3. Further, a compilation by Lasinski et al. [19] is used to supplement our 
own compilation for 7r ± — p scattering at low momentum. The compiled /3 values and the 
parametrized momentum dependences are shown in Figure [9] for p — p and n — p scattering, 



and Figure 10 for — p scattering. 

The resulting total cross sections for nucleon-nucleus scattering (beryllium and aluminum) 



are shown in Figure 11 The calculated values of otot are compared with measurements of 
Otot for n-Be data. The model predictions agree with the data to within several percent, 
and indicate that (Jtot for proton-nucleus and neutron-nucleus interactions are very similar 
except at the lowest energies, as expected from isospin symmetry. The success of the model 
in reproducing otot neutron-nucleus is taken as an indication that the model can be used for 
proton-nucleus and pion-nucleus interactions, where such a check with data is not possible. 
The spread in values between the data and the model is considered a source of systematic 
uncertainty. 
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FIG. 12: Total hadronic cross sections for -n^-Be calculated using the Glauber model (points) 
for 7T + (left) and 7T~ (right). The Breit-Wigner parametrization based on the Carroll data [37] on 
the A(1232) resonance is shown as a dotted line, while the parametrization of the Glauber model 
points is shown as a solid line. The Geant4 default model is shown as a dashed line. 
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FIG. 13: Total hadronic cross sections for tt^-AI calculated using the Glauber model (points) for 
7T + (left) and 7r~ (right). The Breit-Wigner parametrization based on the Carroll data [37] on 
the A(1232) resonance is shown as a dotted line, while the parametrization of the Glauber model 
points is shown as a solid line. The Geant4 default model is shown as a dashed line. 



The <Jtot values obtained for pion-nucleus interactions are shown in Figure 12 for 7T ± -Be 
interactions and Figure 13 for 7r ± -Al interactions. The calculated points are parametrized 
by the black curve. At low momentum (< 600MeV/c for beryllium, < 800MeV/c for alu- 
minum), where the A resonance dominates the cross section, parametrizations based on 
gtot measurements by Carroll et al. are used [H7]. While not used in the flux prediction, 
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FIG. 14: p-Be (left) and p-Al (right) inelastic cross sections measured from Gachurin et al. [38] 
(1 — 4 GeV/c) and Bobchenko et a/. [39 (5 — 9 GeV/c). The solid line is the parametrization used in 
the flux prediction, while the dashed line shows the Geant4 default parametrization. 



the otot values used in the GHEISHA model (the Geant4 default) are shown as a dashed 
black line for comparison. 



B. Inelastic and Quasi-elastic Processes 

In the case of inelastic scattering {ojne\ a much larger set of cross section measurements 
exists eliminating the need for theoretical models. The entire momentum range of interest 
for nucleon-nucleus inelastic scattering and a large subset of the momentum range for pion- 
nucleus inelastic scattering has been measured. 

The available measurements of <Jine for p-Be and p-Al interactions in the relevant mo- 



mentum range are shown in Figure 14 The Gachurin et al. data [38J spans the low momen- 
tum region, while the Bobchenko et al. [39j data covers the high momentum region up to 
9 GeV/ c. Together, they cover the entire momentum range of interest for MiniBooNE. The 
parameterization used to model the momentum dependence is shown as a solid black line. 



Likewise, (Tine for ^-nucleus interactions is shown in Figure 15 The Ashery et al. [4*0] 
measurements are used around the A resonance, while the Allardyce et al. [JT], Gachurin et 
al. [38] , and Bobchenko et al. [39j data are used at higher momentum. The low momentum 
data do not include beryllium; for these points, the cross sections are extrapolated using 
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FIG. 15: Inelastic cross sections for 7r + -Be (top left), 7r~-Be (top right), 7r + -Al (bottom left) and 
7r~-Al (bottom right) as measured in References [40^ (squares), |41J (triangles) and [39 j (circles). 
The solid lines are the parametrizations used in the flux prediction, while the dashed lines are the 
default Geant4 parameterizations. 



the cross sections measured on different elements at the same momentum. The measured 
cross sections are parametrized as A n based on the A-dependence of the cross sections, with 
typical values of n ranging from 0.6 to 0.8. The resulting function is used to infer the cross 
section at A — 9. 

A subset of the inelastic interactions results from quasi-elastic scattering, where hadrons 
scatter off the individual nucleons in the nucleus in a manner analogous to hadron elastic 
scattering off free nucleons. The rate of this process relative to other forms of inelastic scat- 
tering is important since it allows the incoming hadron to emerge from inelastic scattering 
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FIG. 16: Calculated values of ctqel for p-Be (left) and p-Al (right) interactions along with the 
parametrization used in the flux prediction. 



with its initial momentum largely intact, whereas it would otherwise have been effectively 
absorbed or significantly reduced in momentum. 

Unfortunately, the available measurements of ctqel are sparse, with only a few measure- 
ments for pions at low momentum. As a result, we must appeal to a theoretical calculation 
for this part of the inelastic cross section. This can be effected via the shadowed multiple 
scattering expansion, in which oqel is calculated as the incoherent sum of the cross section 
for hadrons to scatter off the nucleons in the nucleus, accounting for the attenuation of the 
hadron wavefunction as it traverses through the nucleus [23] . The cross section for multiple 
scattering of the hadron within the same nucleus can also be calculated in this formalism. 
This is found to be a small fraction of the single-scatter cross section in our case. 

The calculated values of oqel for nucleon-nucleus quasi-elastic scattering are shown in 



Figure 16, while the values for pion-nucleus scattering from Reference [30] are shown in 
Figure 17 The latter figure includes measurements of <Jqel for ^-nucleus interactions 
around the A resonance. The calculated values, along with the measurements, have been 
incorporated into the parametrizations of the momentum dependence of oqel for each of the 
hadron/nucleus combinations. As before, these measurements do not exist for beryllium and 
have been extrapolated assuming an A n dependence, where n has been determined from the 
A-dependence of the measured cross section at each momentum from Reference [30]. The 
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FIG. 17: Quasi-elastic cross sections for 7r + -Be (top left), 7r~-Be (top right), 7r + -Al (bottom left) 
and 7r~-Al (bottom right) as measured in References [40 j (squares) and calculated using the shad- 
owed scattering model (circles). The solid black lines are the parametrizations used in the flux 
prediction. 

resulting values of n range from 0.3 to 1.0, varying with incident pion momentum. 
C. Explicit forms for the Cross Section Parametrizations: 

In summarizing the momentum dependence of the nucleon and pion cross sections, we 
have made use of the following form: 



a = A + B x p n + C x In 2 p + D x lnp 



(9) 
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A 


B 


n 


C 


D 


(p/n)-Be 

V"/ / 












(TTOT 


307.8 


0.897 


0.003 


-2.598 


-4.973 


&INE 


186.7 


104.3 


-1.039 


10.38 


-15.83 


&QEL 


164.8 


-40.09 


0.408 


21.40 


-61.45 


(p/n)-Al 












otot 


760.3 


-0.056 


2.485 


6.173 


-41.60 


cine 


470.9 


-0.259 


2.429 


48.86 


-87.19 


&QEL 


255.7 


8.792 


0.0024 


32.24 


-155.9 



TABLE III: Parameter values for proton and neutron on beryllium and aluminum cross sections 
using Equation [9] 



where p is the momentum of the incident particle in GeV/ c. While this form is inspired by 
Regge theory [12]) it is used as a purely empirical description of the cross section. No physical 
significance is attributed to the parameters apart from the ability of the parametrization to 
describe the measured or calculated cross sections with the appropriate parameters. The 
parameters used in the flux prediction for proton and neutron hadronic cross sections on 
beryllium and aluminum are given in Table |III| 

For (Jine and (Jqel in pion-beryllium/aluminum scattering, a more complicated form is 
needed in order to describe the peak in the cross section near the A resonances: 

-m(p)T R 



a = N, 



R 



(10) 



M| — m(p) 2 + im{p)T r 
+ [ 1 + taoh(0,(p - O )) ] x (A + Bp n + C\n 2 p) 
where p is the momentum of the incident particle in GeV/c. The first term describes a 
Breit-Wigner resonance, where m(p) is the invariant mass of the pion/target nucleon system 
in GeV/c 2 assuming that the target is a nucleon at rest. The second term is a simplified 
version of Equation [9] with a threshold behavior described by a hyperbolic tangent function. 
The threshold function allows the second term to dominate at pion momenta above the 
A(1232) resonance. Here also, the approach is purely empirical; the parameters, including 
the resonance terms are extracted in such a way to reproduce as closely as possible the 
measurements, without assigning any physical significance to any of the parameters. In 
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.4 
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Nr Mr Tr 


7r+-Be 














cttot 


0.814 3.418 


237.6 


111.3 


-4.186 


-9.792 




<?ine 


0.400 5.142 


162.3 


-99.79 


-2.407 


-0.423 850.3 1.201 0.375 


®qel 


0.635 3.784 


-2.38 


-81.84 


-2.702 


3.173 


379.9 1.201 0.558 


7T+-A1 














otot 


0.931 3.186 


569.1 


511.3 


-3.79 


-18.50 




cine 


0.295 2.307 1537.4 


-1109.4 


0.057 


14.40 


510.7 1.189 0.185 


&QEL 


0.698 2.134 


40.38 


89.20 


-1.575 


0.335 


229.4 1.189 0.187 





#0 6s 


A 


B 


n 


C 


Nr Mr Tr 


7r -Be 














otot 


0.814 3.418 


237.6 


111.3 


-4.186 


-9.792 




pine 


0.600 2.874 


92.66 


112.2 


-0.486 


7.500 


371.5 1.201 0.233 


&QEL 


0.626 2.504 


-1.559 


46.41 


-0.633 


1.874 


189.0 1.201 0.185 


7T--A1 














otot 


0.931 3.186 


569.1 


511.3 


-3.79 


-18.50 




vine 


0.706 1.685 


997.8 


-457.8 


0.611 


233.4 


446.8 1.189 0.305 


CQEL 


0.633 2.199 


32.52 


85.15 


-1.225 


1.383 


129.1 1.189 0.305 



TABLE IV: Parameter values for ^-(Be/Al) hadronic cross sections. For cjine an d ctqel, 



Equation 10 is used. For gtot, the parametrization of Carroll et al. [37] is used at low momentum, 
while Equation [9] with a threshold term is used at high momentum. 



particular, the various A resonances are not modeled individually. The parameters used in 



the flux simulation using Equation 10 are given in Table |IV| 

As mentioned in Section IV A , the parameterization of Carroll et al. (37] is used for the to- 
tal cross sections on 7r =1= scattering on beryllium and aluminum for momenta up to 600 MeV/ c 
in the former case and 800 MeV/ c in the latter. At higher momentum, the second term of 



Equation 10 is used with the parameters shown in Table IV 
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V. SECONDARY PARTICLE PRODUCTION CROSS SECTIONS 



The primary source of the neutrino flux at MiniBooNE is the decay of secondary particles 
produced in p-Be interactions. The knowledge of the neutrino flux thus critically depends 
on the understanding of the meson production in p-Be interactions. Most of the flux at 
the MiniBooNE detector comes from tt + — > yU + + decays, while the v e flux is dominated 
by three-body decays of kaons (K + and K®) and the decay of muons (primarily produced in 
the decay of pions). The tables in the Monte Carlo simulation describing the double differ- 
ential cross sections which specify the multiplicity and kinematic properties of the protons, 
neutrons, 7r ± , K ± and K® produced in p-Be interactions at 8.89 GeV/c are based on hadron 
production measurements with similar kinematic configurations wherever possible. In the 
case of 7r =1= , K + and K° production, the double differential cross sections are summarized as 
parametrizations. The parametrizations are evaluated at each point within the table to de- 
termine the corresponding cross sections. For protons, neutrons, and K~ , the cross sections 
are based on a model of hadronic interactions. 



A. Pion Production Measurements 

The cross section tables for 7r =l= production in p-Be interactions are based on parameteriza- 
tions of measurements taken by the HARP [33] and BNL E910 [H] experiments. The HARP 
experiment measured the 7r =l= differential production cross section for p-Be interactions using 
replicas of the MiniBooNE beryllium target at the same incident proton momentum of 8.89 
GeV/c. However, since the analysis of the data from the replica targets is not complete, 
the data used in modeling the pion production is from the thin target run, where a 5% 
interaction length beryllium target was measured. The pion tracks are binned in total pion 
momentum p w ranging from 0.75 to 6.5 GeV/c and angle 9 n with respect to the incident 
proton direction from 30 to 210 mrad. The measurements from the experiment represent 
the average differential cross section over the bin. A complete covariance matrix is also re- 
ported to account for bin-to-bin correlations in the uncertainties. The quoted normalization 
uncertainty is a harp = 2%. 



The left plot of Figure 18 shows the kinematic distribution (in terms of 9^ and p w ) from 



the Monte Carlo simulation of pions produced in the target that decay to produce neutrinos 
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FIG. 18: Left: Production angle 9 vs. momentum p for the tt + in the flux simulation that contribute 
to the Vn flux at the MiniBooNE detector. The color scale indicates the relative cross section for 
7r + production in each bin of angle and momentum. The black box marks the kinematic range 
covered by the HARP measurements [33]. Right: Transverse momentum pt vs. the Feynman- 
scaling variable xf for the K + in the flux simulation that contribute to the neutrino flux at the 
MiniBooNE detector (black squares). The colored points indicate the kinematic regions measured 
by p-Be K + production measurements [3J. 



in the MiniBooNE detector. The black square indicates the kinematic range covered by the 
HARP measurements. 

The BNL E910 experiment measured the tt^ 1 differential cross section for p-Be interac- 
tions at three different energies of the incident protons (6.4, 12.3, 17.5 and 17.6 GeV/c). 
Data is binned in pion momentum p n ranging from 0.4 to 5.6 GeV/c and angle 9^ from 
18 to 400 mrad. The extended coverage of the E910 measurements in the forward angular 
region provides further constraints of the pion production in this kinematic region. A covari- 
ance matrix was not reported for these measurements, hence we use a diagonal bin-to-bin 
covariance matrix. The quoted normalization error is (Jesio — 5%. 

The 7r + production cross sections (momentum distribution of pion production in bins of 



production angle) from the two experiments are shown in Figures 19 (HARP 8.89GeV/c) 



20 (E910 6.4GeV/c), and 21 (E910 12.3GeV/c). The corresponding it production measure- 



ments are shown in Figures 22 , 23 and [24l While a significant body of historical p-Be pion 
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production data exists ([351 H6j 133 133 132]), the measurements are removed from the primary 
beam momentum in the BNB, have insufficient kinematic coverage, or have inconsistencies 
that led to the exclusive use of the latest data from E910 and HARP. The E910 17.5 GeV/c 
and 17.6GeV/c data is also not used for the first reason. 



B. Sanford-Wang Fit to the Pion Production Data 

Following the K2K experiment [50] , the parametrization of Sanford and Wang (SW) 
[5T] is used to describe the 7r ± differential production cross section across different incident 
primary beam momenta. The SW parameterization for the production cross section of a 
given meson species is given by 

— °—{v, 0) = c lP C2 (l exp ( - c; P — 5 - c 6 9(p - c 7 p B cos c « 9)) (11) 

dpdu \ Pb — eg/ \ p^ J 

where i s the double differential cross section, p is the total momentum of the meson in 
GeV/ c, 9 is the angle of the meson with respect to the incident proton in radians, ps is the 
momentum of the incident proton in GeV/ c, and c\ , . . , eg are parameters to be determined 
in the fit to the production data. For the fits to the pion production data, eg is set to unity; 
for the kaon production fits (see below), it is a free parameter. The parametrization allows 
one to relate production data at different incident proton energies, to smoothly interpolate 
the behavior of the cross section between measured points, and to extrapolate into regions 
where production data do not exist. Due to the strong correlation between the C3, C4, and 
C5 parameters, the value of C3 is fixed to unity for the n + production fit. For 7r~, the 
C3 parameter is initially floating, but then fixed to its initial best-fit value when the fit 
is iterated. The correlation results from the limited range of proton momentum covered 
by the measurements (6.4, 8.9 and 12.3 GeV/c). As a result, the data has limited ability 
to constrain the cross section dependence on the proton momentum. The predicted pion 
production properties, however, are not affected by this indeterminacy. 

The values of the parameters are determined from a fit to the tt^ 1 production cross 
section data by minimizing the following x 2 function 



* 2 = £ 



J2 (Di, k - N k T t ) V-i (D jjk - N k Tj) + 



(N k - 1) 



1,3 



-I 



(12) 
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where is the z-th data point for the fc-th data set, T{ is the value of the SW function 
for the kinematic parameters for that data point, Vy * is the bin-to-bin covariance matrix 
for the k-th data set, iV& is relative normalization fit parameter for data set k and is the 
quoted normalization uncertainty for data set k. There are three data sets used in the fit, 
namely HARP and E910 6.4 and 12.3GeV/c. The normalization uncertainties for data sets 
at different beam momenta from the same experiment are the same and varied in common 
across these data sets during the fit procedure. 

Since the data represent the average differential cross section over a range of angle and 
momentum ("bin"), bin-centering corrections must be applied. The bin-centering corrections 
are model-dependent since one must assume how the production cross section varies within 
a bin. Here, the SW parameterization is used to evaluate the bin-centering correction 

where (p^, 9j) is the center of bin (i, j) in the (p, 9) space, SW(pf, Of} is the double differential 
cross section returned by EquationJTTJ and Api and A cos Oj are the bin widths. The MINUIT 
fits are iterated with bin-centering corrections until convergence is achieved. Since we are 
concerned for the most part with pion production at 8.89GeV/c proton momentum, the 
dependence on the proton momentum is not important in predicting neutrino fluxes. 

For the fit to the ir + data, the minimized x 2 /degree of freedom (DOF) using the reported 
experimental uncertainties is 1.8. To obtain parameter uncertainties the fit is performed 
with the covariance matrices scaled by this factor to obtain the parameter uncertainties, 
resulting in an effective x 2 °f unity. The fit parameters are shown in the first row of Table 
|V} The error matrix, shown below the parameters in Table [V] is obtained by varying the 
parameters in such a way that the resulting variations in the SW function cover the spread 
in the data points. This corresponds to an envelope of parameter variations in which the 
resulting x 2 is within 8.14 of the minimum determined by the fit. While this corresponds to 
a 68% confidence level parameter envelope for 7 parameters, the x 2 difference is set by the 
desire to have the variations cover the deviations of the data points and their uncertainties. 
The normalization factors obtained from the two data sets also are compatible within the 
systematic uncertainties quoted by the two experiments (Nharp = 0.966, Nemo = 1.048). 

Likewise, the fit to the HARP [52] and E910 n~ production data[44j with nominal errors 
resulted in a best-fit \ 2 °f 1.16/DOF. The experimental uncertainties are scaled by this 
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p n +(GeV/c) p^GeV/c) 

FIG. 19: Comparison of HARP tt + production cross section data [13] (circles) versus in bins of 9 n 
from 8.89 GeV/c p-Be interactions and best fit SW model (solid lines). The dashed lines represent 
the uncertainty band resulting from varying the parameters within their correlated uncertainties, 
as described in the text. 



factor to achieve a x 2 /DOF of unity. The resulting parameters and covariance matrix are 



shown in Table |VI[ where the covariances are shown in the upper right triangle of the matrix 
(including the diagnonal terms) and the correlation coefficients are shown in the lower left 
triangle of the matrix. The parametrizations using these best-fit parameters, along with the 
expected variation due to the parameter uncertainties, are shown along with the production 



data in Figures 19 24 



C. K + Production Measurements 

For charged kaons, whose decays result in a significant contribution to the flux at 
high energies as well as the u e flux through the K e3 decay mode, there are no measurements 
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Cl 


C2 


C3 


C4 


C5 


C6 


C7 


C8 


Value 


220.7 


1.080 1.000 


1.978 


1.32 


5.572 


0.0868 


9.686 


Cl 


1707.2 


1.146 




-17.646 


-15.968 


-8.81 


-0.7347 


-60.816 


C2 


0.139 0.0396 




-0.1072 


-0.0993 


0.0325 


0.0007 


-0.0777 


C3 
C4 


-0.554 


-0.699 




0.5945 


0.5049 


0.0655 


0.0025 


0.198 


C5 


-0.582 


-0.751 




0.986 


0.4411 


0.0568 


0.0025 


0.2271 


C6 


-0.469 


0.359 




0.187 


0.188 


0.2066 


0.0047 


0.1031 


C7 


-0.795 


0.157 




0.145 


0.168 


0.462 


0.0005 


0.0641 


C8 


-0.368 


-0.098 




0.064 


0.085 


0.057 


0.716 16.0189 



TABLE V: Extracted Sanford-Wang parameters ci_8 (first row), covariance matrix (upper right tri- 
angle including diagonals terms) , and correlation coefficients (lower left triangle) , for ir + secondary- 
production in p-Be interactions. There are no entries in the covariance matrix for parameter C3, 
which is fixed in the fit due to its large correlation with C5. 





Cl c 2 


C3 


c 4 


C5 


C6 


C7 


C8 


Value 


213.7 0.9379 5.454 


1.210 


1.284 


4.781 0.07338 


8.329 


Cl 


3688.9 7.61 




-15.666 


-17.48 


-11.329 


-0.9925 


-91.4 


C2 


0.636 0.0388 




-0.0437 


-0.0509 


0.0102 


-0.0009 


-0.1957 


C3 

c 4 


-0.889 -0.765 




0.0841 


0.0895 


0.0301 


0.0029 


0.2588 


C5 


-0.917 -0.823 




0.983 


0.0986 


0.0375 


0.0033 


0.3141 


C6 


-0.467 0.130 




0.260 


0.299 


0.1595 


0.0051 


0.1933 


C7 


-0.731 -0.204 




0.447 


0.470 


0.571 


0.0005 


0.064 


C8 


-0.362 -0.239 




0.215 


0.241 


0.117 


0.689 


17.242 



TABLE VI: Extracted Sanford-Wang parameters ci_§ (first row), the covariance matrix (upper 
right triangle including diagonal terms), and correlation coefficients (lower left triangle) for 7r~ 
secondary production in p-Be interactions. There are no entries in the covariance matrix for 
parameter C3, which is fixed in the fit due to its large correlation with C5. 
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p.*, GeV/c 



1 2 3 

Pjc+ , GeV/c 



FIG. 20: Comparison of E910 ir + production cross section data |44j (circles) versus p n in bins of n 
from 6.4 GeV/c p-Be interactions and best fit SW model (solid lines). The dashed lines represent 
the uncertainty band for 68% confidence level for 7 fit parameters. 



from the HARP or BNL E910. As a result, measurements reported by other experiments 
measuring K + production in p-Be interactions at primary beam momenta close to 8.89 GeV/c 



are used |53j [Ml ESI [56j [57J |58, 59j. The measurements are summarized in Table VII 



Since no measurements of K + exist at the 8.89GeV/c BNB primary momentum, we 
employ the Feynman scaling hypothesis to relate K + production measurements at different 
primary energies to the expected production at 8.89GeV/c. Theoretically, Feynman scaling 
should be a better model for comparing data from different primary beam momenta. This 
is born out by comparisons of data scaled to the BNB momentum of 8.89GeV/c as shown 



in Figure [25] The hypothesis states that the invariant cross section is a function of only two 
variables, namely Xp and p?, where 

pCm 

X F = max,cm (14) 

P\\ 

is the Feynman scaling variable, defined as the ratio of the parallel component of the mo- 
mentum of the produced particle in the center-of-mass frame and the maximum possible 
value of this quantity for the given reaction, and pt is the transverse component of the 
momentum of the produced particle. In calculating xp, the p™ 8 ^" 1 value is taken from 
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FIG. 21: Comparison of E910 tt + production cross section data |44| (circles) versus p n in bins of n 
from 12.3 GeV/c p-Be interactions and best-fit SW model (solid lines). The dashed lines represent 
the uncertainty band resulting from varying the parameters within their correlated uncertainties, 
as described in the text. 



the exclusive channel p + (p/n) — > A + (p/n) + K + . A more complete description of the 
Feynman scaling fit procedure and results can be found in Reference [50] • 

The Feynman scaling is used as a basis for parametrizing the production data using the 
variables xp and pr- This motivates a six-parameter model given by: 

d 2 a _ P 2 K + ( P Al\ _ P 2 K + 



exp [-c 2 pr - c 3 \x F \ C4 - c 5 p% - c 7 \p T x x F \ Cfi ] 

The model is basically a translation of the Feynman scaling hypothesis where the invariant 
cross section is only a function of Xp and px- It incorporates an exponentially falling p F 
distribution, correlations between p F and Xp, a flat rapidity plateau at xp = and zero cross 
section as xp — > 1. The kinematic threshold constraint is imposed by setting the function 



equal to zero for \xp\ > 1. Figure 25 shows the momentum distribution of the data scaled 
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FIG. 22: Comparison of HARP tt production cross section data |52] (circles) versus p n in bins of n 
from 8.89 GeV/c p-Be interactions and best-fit SW model (solid lines). The dashed lines represent 
the uncertainty band resulting from varying the parameters within their correlated uncertainties, 
as described in the text. 



to 8.89 GeV/c primary momentum in bins of scaled K + production angle. The data include 
normalizations factors obtained from the fit procedure described below. The right plot of 



Figure 18 shows as boxes the xp versus px distribution of K + produced in the target that 
produce neutrinos at the MiniBooNE detector. The colored points indicate the kinematic 
coverage of the various measurements in these two variables. 

The q parameters are determined in a x 2 fit to the production data for 1.2 < pfJ^ 3 < 5.5, 
where p^+ B is the kaon momentum translated to the BNB primary energy using Feynman 
scaling. The p^+ B requirement eliminates most of the data at negative xp, where nuclear 



effects are expected to be dominant. The x 2 takes the same form as in Equation 12, where 
the covariance matrix from the the experiments is diagonal, and the quoted normalization 
uncertainties are used to constrain the normalization factors. The Vorontsov data [59J has 
indications of an error in the normalization outside of their quoted uncertainties. As a 
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FIG. 23: Comparison of E910 7r _ production cross section data |44j (circles) versus p n in bins of 6 n 
from 6.4 GeV/c p-Be interactions and best-fit SW model (solid lines). The dashed lines represent 
the uncertainty band resulting from varying the parameters within their correlated uncertainties, 
as described in the text. 



result, a large normalization uncertainty (500%) was assigned to these measurements with 
the effect that the measurements from this experiment contributes only "shape" information 

since 



without any normalization constraint. The discrepancy is not apparent in Figure 25 



the fitted normalization factor is applied to the measured cross sections. The Lundy data 
[57] were excluded from the fit due to inconsistencies with the other measurements, while 
the Marmer measurements [58] were excluded by the p^+ B requirement. 

The x 2 /DOF f or the fit is 2.28. The errors in the measurements are inflated by a factor 



of = 1.5 to bring the x 2 /DOF to unity. The covariance matrix of the parameters is 

summarizes 



extracted in the same way as employed for the pion production fits. Table VIII 



the results, with the first row listing the seven best-fit parameters, and the 7x7 matrix 
below it listing the covariance. While the parametrization represents the measurements 



quite well (as shown in Figure 25), the uncertainties are inflated by a further factor of four 
for the z/ M — > v e oscillation analysis to account for inconsistencies within the production 
data, the use of the Feynman-scaling hypothesis to relate the production measurements 
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FIG. 24: Comparison of E910 ir production cross section data |44| (circles) versus p n in bins of n 
from 12.3 GeV/c p-Be interactions and best-fit SW model (solid line). The dashed lines represent 
the uncertainty band resulting from varying the parameters within their correlated uncertainties. 



from experiments with proton momenta different from the 8.89GeV/c used at the BNB, 
and a possible discrepancy in the rate of events observed in the MiniBooNE detector 
compared with the predictions based on the beam line simulation|3|. The covariance matrix 



in Table |VIII| and the uncertainty bands in Figure [25] do not include this further inflation of 
the uncertainties by a factor of four. 



Figure 26 shows a similar comparison of the E910 and HARP 7r + production data where 
the E910 data have been scaled to 8.89 GeV/c primary momentum. The data indicates that 
the 7r + production is also consistent with the Feynman scaling hypothesis. When fit to 



Equation 15, a slightly poorer x 2 results than in the SW fit. As a result, the Feynman 
scaling model is not used in the 7r + production model. 
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Dataset 


PhcCLTTL 


Pk+ 


e K + 


Xp 


Pt Civ 




( GeV/c) 


( GeV/c) 


(degrees) 




( GeV/c) 


Abbott 153] 


14.6 


2-8 


20-30 


-0.12-0.07 


0.2-0.7 10% 


Aleshin 09] 

t i 


9.5 


3-6.5 


3.5 


0.3-0.8 


0.2-0.4 10% 


Allaby [51] 


19.2 


3-16 


0-7 


0.3-0.9 


0.1-1.0 15% 


Dekkers |55| 


18.8, 23.1 


4-12 


0, 5 


0.1-0.5 


0.0-1.2 20% 


Eichten 56 

1 ! 


24.0 


4-18 


0-6 


0.1-0.8 


0.1-1.2 20% 


Lundy [57] 


13.4 


3-6 


2, 4, 8 


0.1-0.6 


0.1-1.2 20% 


Mariner [58] 


12.3 


0.5-1.0 


0, 5, 10 


-0.3-1.0 


0.15-0.5 20% 


Vorontsov [59] 


10.1 


1-4.5 


3.5 


0.03-0.5 


0.1-0.25 25% 



TABLE VII: Summary of K + production measurements in p-Be interactions used to characterize 
K + production in the BNB. The table includes -Pf, eam , the primary proton momenta in the mea- 
surement, the momentum and angular ranges of the measurements, as well as the corresponding 
ranges of the Feynman scaling variable xp and transverse momentum pp. Finally, the quoted 
overall normalization uncertainty is listed. 

D. Production of K°, K~ and other particles 

A scheme similar to that used to parametrize the 7r + and K + production data is used for 
neutral kaon production, for which the K e % decay mode of the K® is a source of background 
v e . Since the kaons are produced in strong interactions as K° and K° (primarily the former), 
the kaons have equal content as K® and K\. As a result, the production properties of neutral 
kaons decaying as K® can be obtained by measuring the K® production properties. While 
the K® can contribute to the neutrino flux via the decay of the charged pions produced in 
the — ► 7r + + n~ decay, the most important consideration is the production of v e from 
the decay of the K®. The long life time of the K®, together with the fact that they are not 
focussed, lead to the expectation that the contribution of neutrinos for this source will be 
small relative to the K + . 

The primary source of data for the parametrization comes from two measurements of K® 
production in p-Be interactions in the BNL E910 experiment (pbeam = 12.3 and 17.5 GeV/c) 
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FIG. 25: Comparison of K + invariant cross section data (points) as a function of Pk+> the K + 
momentum, in bins of @k+> the K + production angle (in radians), with the Feynman-scaling based 
parametrization with best-fit parameters shown as a solid line. The scaling has been used to relate 
the measurements at different primary beam momenta to the 8.89 GeV/ c primary momentum in the 
BNB. Normalization factors from the best-fit are applied to the data. The dashed lines represent 
the uncertainty band resulting from varying the parameters within their correlated uncertainties. 
The uncertainty bands include the factor 1.5 error inflation to set x 2 /DOF = 1. 

and the measurements of Abe et al. [61] (12.3GeV/c) at KEK. Since the neutral kaons are 
not focused by the magnetic field of the horn, the forward production (< 5°) is particularly 
relevant for predicting the BNB neutrino flux. While the production data from the BNL 
E910 and KEK measurements do not cover this region, the combination of the two data 
sets are sufficient to constrain the production cross section in this forward region via the 
Sanford-Wang parametrization. The extracted parameter values and covariance matrix are 



45 





Cl 


C2 


C3 


C4 


cs 


C6 


C7 


Value 


11.70 


0.88 


4.77 


1.51 


2.21 


2.17 


1.51 


cy 


1.094 


0.0502 


2.99-10" 3 


-0.0332 


-0.0375 


0.125 


0.0743 


C'2 


0.378 


0.0161 


1.39-10" 3 


-1.44-10" 3 


-0.0126 


0.0322 


0.022 


C3 


0.033 


0.127 


7.47-10- 3 


2.06- 10~ 3 


1.93-10~ 3 


0.0135 


- 3.34-10~ 3 


C4 


-0.540 


-0.193 


0.405 


3.46- 10~ 3 


2.03- 10~ 3 


-4.1M0- 3 


- 6.28-10~ 3 


C5 


-0.297 


-0.822 


0.185 


0.286 


0.0146 


-0.0154 


-0.0244 


C6 


0.280 


0.595 


0.366 


-0.164 


-0.299 


0.182 


0.126 


C7 


0.178 


0.435 


-0.097 


-0.268 


-0.506 


0.741 


0.159 



TABLE VIII: Best-fit Feynman scaling model parameters q from a fit to K + production data (first 
row). The covariance matrix for the parameters with uncertainties inflated by a factor of 1.5 to 
set x 2 /°OF = 1 is in the upper right triangle (including diagonal terms) of the table below the 
parameters, but does not include the additional inflation by a factor of four described in the text. 
The lower left triangle shows the correlation coefficients. 

summarized in Table HXl 

For K~ production, the scarcity of production measurements in the relevant kinematic 
regions motivated the use of the MARS hadronic interaction package [H2] to determine 
the absolute double differential cross sections. The cross sections are obtained by simulating 
8.89 GeV/c p-Be interactions on a thin beryllium target and recording the rate and spectrum 
of outgoing K~ . The expected relative contribution of neutrinos of all species from K~ 
decays is expected to be small. Neutrino flux contributions from semileptonic hyperon 
decays {e.g. A, S, etc.), estimated using a FLUKA|63J simulation, are also negligible. 

Secondary protons and neutrons emerging from the p-Be inelastic interactions are sim- 
ulated based on the predictions of the MARS model, with the exception of quasi-elastic 
scattering, in which case the final state proton kinematics are handled by a custom model. 
The production of all other particle species is handled by the default Geant4 hadronic model. 

The properties of the particle production model are summarized in Table |Xj The table 
shows the average multiplicity per p — Be reaction (denned as inelastic interactions excluding 
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FIG. 26: Invariant pion production cross section from HARP and E910 versus p n in bins of 9- 
The E910 measurements are rescaled to pb = 8.89 GeV/c. 



quasi-elastic scattering), along with the mean momentum and production angle. The ir + 
and 7r~ production occur with similar multiplicities, though the former tends to be harder 
and more forward directed. The larger overall multiplicity for the ir~ is due to the extrap- 
olation of the cross sections to large angles that are not covered in the HARP and E910 
measurements. Since the contribution to the neutrino flux from such pions is small, the 
impact of uncertainty in this extrapolation is suppressed. The kaon production is an order 
of magnitude smaller than the pion production, with K~ production particularly suppressed 
relative to K + and K° production. 



VI. PREDICTION OF NEUTRINO FLUX AT MINIBOONE 



The results of the simulation are summarized in Figures 27 and 28 which show the total 
predicted flux of each neutrino species at the MiniBooNE detector in neutrino mode and 
anti-neutrino modes, respectively. In each case, the z/ e /z/ e contribution is less than 1% at 
the peak of the v^jv^ flux, though it rises at higher energies. As shown, the predicted 
fluxes exhibit many features that are better understood by analyzing the sources of each 
component of the flux. 
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Cl 


C'2 


C.3 


C4 


C5 


C6 


C7 


C8 


eg 


Value 


15.130 


1.975 


4.084 


0.928 


0.731 


4.362 


0.048 


13.300 


1.278 


Cl 


32.3 


-0.09687 


0.8215 


-0.1018 


-0.2124 


-0.8902 


-0.1333 


16.55 


-1.789 


C2 


-0.055 


0.09574 


0.03248 


0.00131 


-0.01303 0.08836 


-0.00031 


-1.536 


-0.2156 


C3 


0.199 


0.144 


0.5283 


-0.01922 


0.02267 


-0.0033 


-0.00236 


0.0391 


-0.08017 


C4 


-0.062 


0.015 


-0.091 


0.08442 


0.00405 


0.00071 


-0.00037 


-0.01443 


-0.07301 


C5 


-0.377 


-0.425 


0.315 


0.141 


0.00982 


0.00287 


0.00028 


-0.05777 


0.02966 


c 6 


-0.261 


0.476 


-0.008 


0.004 


0.048 


0.3599 


0.00385 


-4.751 


-0.1577 


C7 


-0.724 


-0.031 


-0.100 


-0.039 


0.087 


0.198 


0.00105 


0.05806 


0.00686 


C8 


0.255 


-0.435 


0.005 


-0.004 


-0.051 


-0.694 


0.157 


130.2 


1.222 


eg 


-0.183 


-0.406 


-0.064 


-0.146 


0.174 


-0.153 


0.123 


0.062 


2.948 



TABLE IX: Best-fit Sanford-Wang model parameters Cj from a fit to production data (first 
row). The covariance matrix for the parameters is in the upper right triangle (including diagonal 
terms) of the table below the parameters, while the bottom left corner of the table reports the 
corresponding correlation coefficients. 



The integrated contribution of each (anti-)neutrino species, along with their dominant 



decay chains, are shown in Table XI for the neutrino-mode horn configuration, and Table 



XII for the anti-neutrino-mode horn configuration. The dominant contribution from decay 



chains in which the parent meson is produced by a nucleon is separated from those in 
which it is produced by a meson interaction. This is due to the qualitatively different 
level of systematic understanding for the two processes. For the former, the production 
cross sections are based on the particle production experiments described in Section |Vj with 
systematic uncertainties propagated from the uncertainties reported by these experiments. 
For the latter, the simulation relies on the default Geant4 hadronic interaction model to 
provide the production cross sections. Fortunately, the latter is a small contribution to the 
flux in all cases. 



Figure 29 shows the channels through which the and are produced in neutrino mode. 
For the Vn flux, the tt + — > v„ contribution is dominant for energies less than 2 GeV, while 
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FIG. 27: Total predicted flux at the MiniBooNE detector by neutrino species with horn in neutrino 
mode. 




FIG. 28: Total predicted flux at the MiniBooNE detector by neutrino species with horn in anti- 
neutrino mode. 
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,-9 _ v„ channels 



— all 




E v (GeV) 

FIG. 29: Predicted (top) and z7„ (bottom) fluxes at the MiniBooNE detector by parent meson 
species with horn in neutrino mode. The black line is the total predicted flux, while all the 
subcomponents apart from the dashed black are from nucleon-induced meson production of the 
indicated decay chains. The dashed black histogram includes all other contributions, primarily 
from meson decay chains initiated by meson- nucleus interactions. 
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FIG. 30: Predicted u e (top) and V e (bottom) flux at the MiniBooNE detector by parent meson 
species with horn in neutrino mode The black line is the total predicted flux, while all the subcom- 
ponents apart from the dashed black are from nucleon-induced meson production of the indicated 
decay chains. The dashed black histogram includes all other contributions, primarily from meson 
decay chains initiated by meson-nucleus interactions. 
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Particle 


Multiplicity 


(P) 


(0) 




per reaction 


( GeV/c) 


(mrad) 


p 


1.5462 


2.64 


441 


n 


1.3434 


1.59 


586 




0.9004* 


0.82 


556 


7T+ 


0.8825* 


1.11 


412 


K+ 


0.0689 


1.69 


332 


K° 


0.0241 


1.34 


414 


K~ 


0.0024 


1.26 


259 


Total 


4.7679 


1.69 


496 



TABLE X: Average multiplicity per particle-producing reaction for secondary particles produced in 
the inelastic collisions of 8.89 GeV/ c primary protons on beryllium, as well as average momentum (p) 
and angle (9) with respect to the primary proton direction. Multiplicities and average kinematics 
refer to particles produced in the forward hemisphere in the laboratory frame and with transverse 
momentum less than 1 GeV/c. *see comment in text. 



the K + — > flux become dominant at higher energies. The two peaks in the K + flux at 
low energies are from two- and three-body K + decays at rest. Due to the relative size of the 
7r + flux, however, they are not visible in the total flux. There is a small contribution to 
the flux from pions produced in the decay of kaons, and a similar contribution from tertiary 
meson- induced production of other mesons that decay to produce v^. 

For the in neutrino mode, 7r~ — ► flux is dominant at all energies. The next largest 
contribution comes from the 7r + — > fi + — > decay chain. For the flux, the analogous 
contribution from the 7r~ — ► fi~ —>■ decay chain is suppressed by the defocusing of the 
7r~. The kaon contribution is suppressed by the lower rate of K~ production relative to K + 
production. Apart from low energies (< 200 MeV) the predicted flux is typically ~ 6% 
of the flux. 

The channels through which v e and V e are produced in neutrino mode are shown in 

-> v e decay 



Figure |30j For the v e flux, the two dominant components are the n + — ► [i + 
and three-body K + — > v e decay, where the former is dominant at low energies (< 1 GeV) 
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Flux (z//cm 2 /POT) 
Frac. of Total 


5.19 x 10~ 10 
93.6% 


3.26 x 10 -11 
5.86% 


Composition 


7T+: 96.72% 
K+: 2.65% 
K+ -► 7T+: 0.26% 
k° — v 7r+- n n4% 

-TV ^ /I . U.U^/O 

if : 0.03% 
it' ^ 0.01% 
Other: 0.30% 


7T _ : 89.74% 
vr+ -► 4.54% 
A -- : 0.51% 

K°^7T-: 0.24% 
iY+ -► 0.06% 
iY-^vr-: 0.03% 
Other: 4.43% 








Flux (zVcm 2 /POT) 
Frac. of Total 


2.87 x 10~ 12 
0.52% 


3.00 x 10~ 13 
0.05% 


Composition 


it+ -> n+: 51.64% 
K+: 37.28% 
i\T°: 7.39% 
7T+: 2.16% 
->■ n+: 0.69% 
Other: 0.84% 


iY°: 70.65% 
ir~ -► /i - 19.33% 
X-: 4.07% 
7T - : 1.26% 
K~ ^ 0.07% 
Other: 4.62% 



TABLE XI: Predicted v^/Vp (top) and i/ e /i/ e (bottom) fluxes at the MiniBooNE detector with 
horn in neutrino mode. The contribution of flux from meson decays where the parent particle 
in the decay chain is produced by proton or neutron interaction. The "other" category includes 
channels with contributions less than those shown, along with cases where the parent particle in 
the decay chain is produced by a meson interaction. 

and the latter is dominant at higher energies. The peak in the K + — > v e spectrum at low 
energies is from the decay of K + at rest (the peak from two-body decay is much smaller due 
to helicity suppression). For V e , the n" — > fj,~ — > V e flux contributes only at lower energies 
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Flux (zVcm 2 /POT) 
Frac. of Total 


5.42 x 10" 11 
15.71% 


2.93 x 10" 10 
83.73% 


Composition 


vr + : 88.79% 
K+: 7.53% 

7T~ — > n - - 1 77% 

if : 0.26% 
Other: 2.00% 


Tr~: 98.4% 
K~: 0.18% 
R' — > 7T _ - n 0^% 

XV ^ /I . u.uu/o 

K°: 0.05% 
vr+ -> ^+: 0.03% 
K-^7T~: 0.02% 
Other: 1.30% 








Flux (zVcm 2 /POT) 
Frac. of Total 


6.71 x 10~ 13 
0.2% 


1.27 x 10~ 12 
0.4% 


Composition 


K+ : 51.72% 
ET° : 31.56% 
7T+ -> 0+ 13.30% 
7T+ : 0.83% 
-> 0.41% 
Other: 2.17% 


7T _ 75.67% 
K°: 16.51% 
if": 3.08% 
7T _ : 2.58% 
K~ ^fT: 0.06% 
Other 2.10% 



TABLE XII: Predicted v^/v^ (top) anci ^e/^e (bottom) fluxes at the MiniBooNE detector with 
horn in anti-neutrino mode. The contribution of flux from meson decays where the parent particle 
in the decay chain is produced by proton or neutron interaction. The "other" category includes 
channels with contributions less than those shown, along with cases where the parent particle in 
the decay chain is produced by a meson interaction. 

due to the defocusing of the 7r~, and the K~ — > V e contribution is suppressed both by the 
lower production rates and the defocusing. The rest of the spectrum is dominated by K® 
decay. As in the v^jv^ case, the predicted V e flux is ~ 10% of the v e flux. 



Figures 31 and 32 show a similar composition for the predicted anti-neutrino mode flux. 
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FIG. 31: Predicted (top) and z7„ (bottom) fluxes at the MiniBooNE detector by parent meson 
species with horn in anti-neutrino mode. The black line is the total predicted flux, while all 
the subcomponents apart from the dashed black are from nucleon-induced meson production. 
The dashed black histogram includes all other contributions, primarily from meson decay chains 
initiated by meson- nucleus interactions. 
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E v (GeV) 

FIG. 32: Predicted u e (top) and v e (bottom) fluxes at the MiniBooNE detector by parent meson 
species with horn in anti-neutrino mode. The black line is the total predicted flux, while all 
the subcomponents apart from the dashed black are from nucleon-induced meson production. 
The dashed black histogram includes all other contributions, primarily from meson decay chains 
initiated by meson- nucleus interactions. 
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The flux is dominated at all energies by 7r~ — > z/ M decays; the suppressed production 
of K~ results in the contribution of K~ — > being much smaller than the corresponding 
K + — > z/^ contribution in neutrino mode. Furthermore, since the K~ that come to rest are 
captured, the flux does not show the peaks from two-body and three-body decay-at-rest 
at low energies that are found in the from K + decay in both neutrino and anti-neutrino 
mode. It can also be seen that the high energy flux of is not substantially suppressed 
relative to the V^. In fact, despite the defocusing of K + , the K + — > flux is larger than that 
of the K~ — > decay. This is due to the relative production rates and, at high energies, the 
leading particle effect where n + and K + have a harder momentum spectrum relative to their 
negatively-charged counterparts. The high momentum of the particles that produce these 
neutrinos, along with their forward angular distribution, result in less defocusing from the 
horn for the wrong-sign component (positive (negative) particles for (anti-)neutrino mode). 
A similar effect is seen for the v e jv e components in anti-neutrino mode: while the V e are 
dominated by fi~ decays at energies below 2 GeV, the K + — > v e flux is larger than the 
K~ — > T> e flux. A related observation is the fact that while the absolute rate of v e jv e 
from K® is unchanged from neutrino mode, the relative contribution is much stronger in 
anti-neutrino mode. 



VII. SYSTEMATIC UNCERTAINTIES 



The systematic uncertainties in the neutrino flux prediction come from several sources: 

• Proton delivery: The simulation determines the rate and spectrum of neutrinos per 
proton-on-target. This information is combined with the number of protons delivered 
to the target to determine the number of neutrinos passing through the MiniBooNE 
detector over the data collection period. As a result, the predicted number of neu- 
trino interactions in the detector varies directly with the uncertainty in the number 
of protons-on-target. A related uncertainty arises from the optics of the proton beam 
which can change the expected number of protons interacting in the target (or else- 
where), changing the neutrino flux. 

• Particle Production: The uncertainties in the rate and spectrum of secondary particles 
produced in the p-Be interactions likewise affect the rate and spectrum of the neutrinos 
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they produce. This is the dominant uncertainty. 

• Hadronic Interactions: The rate of hadronic interactions affect many aspects of the 
neutrino production, including the rate of p-Be interactions as well as the probability 
for mesons to survive possible hardonic interactions in the target or horn and decay 
to produce neutrinos. Uncertainties in the rate of these interactions affect both the 
rate and shape of the flux. 

• Horn magnetic field: The focusing properties of the horn change with the current as 
well as the distribution of the magnetic field within the conducting elements. Uncer- 
tainties in these properties result in spectral distortions of the neutrino flux. 

• Beamline geometry: Misalignments or displacements of the beamline components from 
their expected orientation and locations can affect the neutrino flux in many ways. For 
example, a misalignment can result in the detector being exposed to a different part 
of the neutrino flux than expected. A displacement of the target with respect to the 
horn can result in a variation in the focusing properties. 



A. Proton delivery 

The systematic uncertainties associated with the delivery of the primary proton beam 
to the beryllium target can be divided into two parts: the uncertainty in the number of 
protons delivered to the beamline and the uncertainty in the number which actually strike 
the target. Having entered the target, there are further uncertainties associated with how 
often the protons will interact to produce secondary particles based on the assumed hadronic 



cross sections; we consider these uncertainties in Section [VII C 

As mentioned in Section [TTJ, the protons delivered to the BNB are measured by two 
toroids upstream of the target. The systematic uncertainty in the resulting spill-by-spill 
measurements has been estimated to be 2% based on uncertainties in the toroid circuit 
elements and uncertainties in the calibration procedure. Since the overall neutrino flux scales 
with the delivered protons, this source or error can be treated as an overall normalization 
uncertainty. The toroid measurements have been cross checked by measuring the activation 
on a gold foil inserted into the beam. The number of protons striking the foil inferred from 
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this measurement agree with the toroid measurements within the ~ 10% uncertainty of the 
measurement. 

The effect of uncertainties in the primary beam optics, most notably the transverse 
profile and focusing and divergence properties, have been estimated by simulating the effects 
introduced by perturbing the default beam parameters. A number of different configurations, 
including varying the focal point across the length of the target, a "pin" configuration with no 
transverse spread or angular divergence, and a "pencil" configuration with transverse spread 
but no angular divergence, have been considered. The resulting changes to the number of 
protons expected to interact in the target is less than 1%, which is taken as a systematic 
uncertainty in the overall normalization of the neutrino flux. 



B. Particle production 

The uncertainties in the particle production are summarized as a covariance matrix in 
the fitted parameters of the functions parametrizing the double differential cross section as 
described in Section |Vj The effect of these uncertainties is propagated to the neutrino flux 
by drawing random parameter vectors according to the covariance matrix via the Cholesky 
decomposition [TS] . The resulting variation in the double differential meson production cross 
section at any point in (p, 6) can be evaluated with respect to the default value. The change 
in the neutrino flux can then be recalculated by assigning a weight corresponding to the 
ratio of the double differential cross section of the secondary particle producing the neutrino 
with the varied and default parameters. 

In this way, the flux resulting from different production distributions summarized by 
alternate parameters can be calculated without re-running the flux simulation. By accumu- 
lating the covariance of the flux distribution as the parameters are varied according to their 
covariance matrix, the uncertainties are propagated into the neutrino flux. This procedure 
is repeated for each parent particle species n~ , K + , K®), and for each neutrino species 
(up, v e , V e ) to obtain the total flux uncertainty, accounting for the correlated variations 
in the different neutrino species. This results in a covariance matrix for the predicted flux 
of each neutrino species from each of the meson species. 



Figure 33 shows the fractional uncertainty in the neutrino flux from n + and K + produc- 



tion uncertainties, corresponding to the square-root of the diagonal entries of the covariance 
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FIG. 33: Left: The fractional uncertainties in the tt + — ► and 7r + — ► /U 



> i/ e flux with the horn 

in neutrino mode due to uncertainties in the ir + production in p-Be interactions. Right: Same for 
the K + — » Vn and K + — > v e flux from uncertainties in the K + production in p-Be interactions. 
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FIG. 34: Left: The fractional uncertainties in the tt~ — > and 7r _ — > n~ — ► V e flux with the horn 
in anti-neutrino mode due to uncertainties in the ir~ production in p-Be interactions. Right: Same 
for the K + — > i/„ and K + — > v e flux from uncertainties in the K + production in p-Be interactions. 
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Energy (GeV) v e Energy (GeV) Energy (GeV) 



FIG. 35: Left: Correlation matrix for variations in the tt + — > v„ flux due to uncertainties in 
7r + production in p-Be interactions with the horn in neutrino mode. Center/Right: Same for 
variations in the tt + — ► /x + — > z/ e flux (center) and correlations between the tt + — ► i/„ flux and the 
7r + — > fi + — > i/ e flux (right). The color scale on each plot ranges from to 1. 




Energy (GeV) v e Energy (GeV) Energy (GeV) 

FIG. 36: Left: Correlation matrix for variations in the K + — ► i/„ flux due to uncertainties in K + 
production in p-Be interactions with the horn in neutrino mode. Center /Right: Same for variations 
in the K + — * v e flux (center) and correlations between the K + — * flux and the K + — * v e flux 
(right). The color scale on each plot ranges from to 1. 

matrix resulting from the procedure described above divided by the predicted flux. In the 
left plot, the solid histogram shows the fractional uncertainty in the flux of neutrinos at the 
MiniBooNE detector from 7r + — > produced in p-Be interactions due to the uncertainties 
in the tt + production. The strong correlation between the energy of the and the energy 
of the 7r + which decayed to produce it results in a large rise in the fractional uncertainty at 
neutrino energies greater than 2 GeV reflecting the large uncertainties in high-momentum 
pion production. Likewise, the uncertainty rises at low neutrino energies (< 200 MeV) due 
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FIG. 37: Left: Correlation matrix for variations in the ir~ —* flux due to uncertainties in n~ 
production in p-Be interactions with the horn in anti-neutrino mode. Center/Right: Same for 
variations in the tt~ — > pT — > v e flux (center) and correlations between the tt~ — ► flux and the 
tt~ — > [i~ — > z7 e flux (right). The color scale on each plot ranges from to 1. 
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FIG. 38: Left: Correlation matrix for variations in the K + — > i/„ flux due to uncertainties in 
K + production in p-Be interactions with the horn in anti-neutrino mode. Center/Right: Same 
for variations in the K + — > v e flux (center) and correlations between the K + — > flux and the 
K + — > i/ e flux (right). The color scale on each plot ranges from to 1. 

to the rise in the uncertainties for low momentum 7r + production. Fortunately, relatively 
few neutrinos are produced in this region by the 7r + decays; in the region below 1 GeV where 
the 7r + — > z/ M contribution is dominant, the uncertainty is approximately 17%. 

The dashed histogram shows the fractional uncertainty in the v e flux from the tt + — > 
fi + — > v decay chain resulting from the uncertainties in the 7r + production, the primary 
channel for low energy (< 1 GeV) v e flux. Since the correlation between the energy of the 
v e and the tt + which produces it is weak due to the three-body decay of the muon, the 
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uncertainties are more uniform as a function of energy. 



The right plot in Figure 33 likewise shows the fractional uncertainty for the flux of and 
v e resulting from the decay of K + produced in p-Be interactions due to the uncertainties in 
the K + production. These channels are the primary contribution for with energy greater 
than 2.3 GeV and v e with energy greater than 1.2 GeV. Due to the larger K + mass, the 
correlation between the momentum of the K + and the neutrinos from its decay is also weak. 



Figure [34] shows the corresponding plots for the horn in anti-neutrino mode, where the 
7r~ — > flux is dominant. While the corresponding charged kaon channel would be K~ — > 
(p^/ve), the K + — > (z^ t /V e ) uncertainties are shown instead, since the contribution of this 
channel is larger. 



Figure 35 shows the bin-to-bin correlations in the uncertainties related to the pion pro- 
duction. The left and center plots show the correlation matrix associated with the fractional 
uncertainties in the ir + — > and 7r + — > fi + — > v e flux in p-Be interactions, respectively. 
The 7r + — > flux exhibit correlations that are strongest between nearby bins, with the 
correlations steadily weakening for bins separated by more than several hundred MeV. The 
7r + — > fi + — > v e flux, however, shows correlation between energies which are more widely 
separated, as would be expected from the three-body decay of the /i + that produces this 



flux. The right plot in Figure 35 shows the correlations between the uncertainties in the 
two components of the 7r + flux. As expected, the ir + —>■ flux at a given energy is most 
strongly correlated with 7r + — > fi + — > v e flux at lower energies. 



Figure 36 shows similar correlations for the K + — > (left) and K + — > u e (center) fluxes. 
The situation is quite different from the ir + flux; the uncertainties are correlated across 
energies for each flux. Likewise, the right plot, which shows the correlations between the 
uncertainties in the two sources of neutrinos, is also strongly correlated. This reflects the 
large normalization uncertainty assigned to the K + production; the variations in the K + 
production correspond mainly to shifts in the overall rate. 



Figures [37] and [38] show the corresponding plots for the horn in anti-neutrino mode. 
Once again, the correlations for the K + — > (z/ M /V e ) uncertainties are shown instead of the 
corresponding K~ —>■ {p^/V^j. The results show a similar pattern of correlations to that 
observed in the neutrino mode. 

Due to the large uncertainties in the flux prediction from particle production, which result 
not from uncertainties in the measured particle production but from the parameterizations 
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Aaxor (mb) 


(mb) 

± IV ±J \ J 


A-vqel (mb) 




Be Al 


Be Al 


Be Al 


(p/ra)-(Be/Al) 


± 15.0 ±25.0 


±5 ±10 


± 20 ±45 


7r ± -(Be/Al) 


± 11.9 ±28.7 


± 10 ±20 


± 11.2 ±25.9 



TABLE XIII: Cross section variations for systematic studies. For each hadron-nucleus cross section 
type, the momentum-dependent cross section is offset by the amount shown. 

used to describe the measurements, an alternative description of the measurements was 
investigated. In this method, the particle production data are interpolated as splines using 
the DCSPLN routine from CERNLIB [13]. The two-dimensional measurements of the double 
differential cross section in p and 9 are interpolated as splines first in 9 at fixed values of p. 
The resulting splines in 9 are then interpolated to produce values as a function of p. The 
splines also extrapolate the double differential cross section into region where measurements 
do not exist, most notably at low pion momentum (< 700MeV/c). 

Uncertainties in the particle production are derived by varying the measured double dif- 
ferential cross sections according to the 78 x 78 covariance matrix describing the uncertainties 
in the measurements. Each variation results in an alternate set of double differential cross 
sections that is interpolated. By repeating this process (in practice one thousand times) a 
covariance matrix describing the pion production variations can be derived that can be used 
to propagate the uncertainties into the predicted neutrino flux. As a result of this procedure, 
the pion production uncertainties can be derived in a manner that is directly connected with 
the experimental uncertainties, circumventing the difficulties associated with parameteriz- 
ing the double differential cross section and reducing the uncertainties in the neutrino flux 
significantly. In the core of the energy distribution (0.5 — 1.0 GeV) the flux uncertainty 
resulting from this method is 5 — 7%. At low energies (< 500 MeV), where measurements do 
not exist, the uncertainties rise to ~ 20%. At higher energies > 1 GeV, where the measure- 
ments have larger uncertainties, the flux uncertainty rises to ~ 10%. This method was 
not used for the neutrino oscillation results reported in Reference [2] but illustrates how the 
additional uncertainties associated with the parameterization of particle production mea- 
surements can be avoided in future analyses. 



64 




FIG. 39: Observed q 2 distribution from Bellettini et al. [M] for 20 GeV/c p-Be scattering. The solid 
line represents the fit to two exponential distributions with the dashed and dotted lines showing 
the elastic and quasi-elastic contributions, respectively. 



C. Hadronic Interactions 

Uncertainties due to hadronic interactions are considered by varying the components of 
the hadronic cross sections. First, the total hadronic cross section <Jtot, the total inelastic 
(Tine and the quasi-elastic oqel cross sections are separately varied for nucleons on beryllium 
and aluminum. Second, the same is done for the pion cross sections. In each case, the 
variations are a flat, momentum-independent offset. Due to the various relations between 
the cross sections, the variation of <Jtot results in a variation of oela {?ine is fixed). When 
°ine is varied, the balance between oela and (Tine is changed, while keeping their total at 
otot- Finally, when oqel is varied, the relative proportion of (Jqel to the cross section for 
all other inelastic processes is changed, while keeping their sum (ctine) fixed. 

The variations for otot are based on the agreement of the Glauber model calculations 
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with the available n-nucleus measurements, shown in Figure 11 The deviations of the mea- 
surements from the model are used to set the magnitude of the variation. These variations 
are also applied to otot for ^-nucleus measurements by scaling the variations from the 
nuclear case by the ratio of otot at high momentum in the two cases. Since there are no n - 
nucleus otot data above the A (1232) resonance to verify the model, this assumes that the 
model works as well for ^-nucleus interactions as it does for nucleon-nucleus interactions. 

The variations for (Jine are similarly set by the deviations of the measurements from 
the parametrization. Here the deviations are smaller since the parametrizations are derived 
directly from the data; the variations are intended to incorporate the uncertainties in the 
measurements. 

The uncertainties for oqel are set by comparing the calculated cross section to the in- 
ferred oqel from the measured q 2 distribution in p-Be scattering from Reference [64] shown 



in Figure 39 This distribution is fit to the sum of two exponentials corresponding to the 
elastic and quasi-elastic scattering components. The fitted slope of the quasi-elastic compo- 
nent is consistent with free nucleon-nucleon elastic scattering and leads a ratio uqel/cela 
of 0.6. Since the free nucleon cross sections and otot do not change appreciably from 
8.89 GeV/c to 20 GeV/c, the value and uncertainty of oela at 8.89 GeV/c can be used to ob- 
tain oqel = 44 ±9 mb at the same beam momentum. This can be compared to the 34.9 mb 
obtained from the shadowed scattering model. An uncertainty of 20 mb is assigned to oqel 
to account for both the difference and the uncertainty from the Bellettini measurement. 
The uncertainty for nucleon-aluminum scattering is obtained by scaling this uncertainty by 
the ratio of the predicted oqel in aluminum and beryllium. The variation for oqel in 
nucleus scattering is obtained by scaling the variation in nucleon-nucleus scattering by the 
ratio of oqel for 7r =l= and nucleons. 

Table IXIIII summarizes the variations in all six hadron-nucleus combinations. Of these 
variations, the variations in oqel have the largest effect. 



D. Horn Magnetic Field Modeling 

Uncertainties on two properties of the horn magnetic field result in systematic uncertain- 
ties in the neutrino flux. The first is the horn current. The commercial current transformers 
(Stangenes Industries 3-0.002) have a rated accuracy of 0.5%. The effect of a 1 kA variation 
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in the nominal horn current (174 ± 1 kA) is simulated to set the systematic uncertainty. 

A second source of uncertainty arises from the modeling of the current within the inner 
cylinder due to the so-called "skin effect" . The skin effect allows for the magnetic field to 
penetrate into the conductor, increasing the effective magnetic field experienced by particles 
traversing through the inner conductor into the horn cavity. This is important for particles 
produced at small angles (particularly high momentum) that barely penetrate into the horn 
cavity before exiting the front of the horn into the collimator region. For these particles, 
the bulk of the magnetic field seen by the particle in this trajectory may come from the field 
within the inner conductor. 

The expected current distribution in a cylindrically symmetric configuration was nu- 
merically evaluated and found to be well-approximated by a current density exponentially 
decreasing from the outer surface of the inner conductor with a decay length set by the skin 
depth (1.4 mm). The magnetic field configuration corresponding to this current density is 
simulated and taken as the default configuration. The simulation is also run without the skin 
effect, simulating the situation where the current density lies entirely on the outer surface of 
the inner conductor of the horn, resulting in no magnetic field penetration. The difference 
between these two configurations, a few percent in the predicted neutrino flux for with 
energies < 1 GeV and up to 18% for with energies ~ 2 GeV, is considered a systematic 
uncertainty. 

E. Geometry Uncertainties 

Variations in the geometric configuration of the beamline are simulated to investigate their 
effect on the neutrino flux. These variations include moving the target position relative to 
the rest of the beamline (in particular the horn), varying the radius of the decay pipe, and 
moving the collimator along the beam axis and changing its aperture. The magnitudes of 
the geometric perturbations which are required to effect a substantial (> 1%) change in the 
flux are well outside of what are considered the tolerances and precision of the constructed 
beamline. As a result, no significant systematic uncertainty is assigned to the beamline 
geometry. 
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TABLE XIV: Variations in the total flux of each neutrino species in neutrino mode due to the 
systematic uncertainties. 
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TABLE XV: Variations in the total flux of each neutrino species in anti-neutrino mode due to the 
systematic uncertainties. 



F. Summary of Systematic Uncertainties 



Tables XIV and XV summarize the variations in the total flux for each neutrino species 
resulting from the systematic uncertainties discussed in this Section. By far the largest 
uncertainty arises from the particle production uncertainties. Much of this uncertainty 
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FIG. 40: Top: Change in the flux from ir + due to the dominant sources of systematic uncertainty 
apart from particle production. Left: Ratio of flux from increased horn current (circles), increased 
nucleon-nucleus quasi-elastic cross section (squares) and increased pion-nucleus quasi-elastic cross 
section (upward triangles) to the default flux. Right: Ratio of flux from decreased horn current 
(circles), decreased nucleon-nucleus quasi-elastic cross section (squares), decreased pion-nucleus 
quasi-elastic cross section (upward triangles), and turning off the skin effect (downward trianges), 
to the default flux. Bottom: Same for predicted flux from ir~ in antineutrino mode. 



arises not from the accuracy of the measurements, but from the parametrizations used to 
summarize the double differential cross sections. These latter contributions manifest as 
inflated uncertainties on the parameters resulting from the x 2 /DOF at the best-fit values 
and other considerations such as the dependence on the choice of parameterization. 
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The flux-averaged uncertainties provide a rough gauge to the relative size of the various 



uncertainties; they are not used in the z/„ — ► v e oscillation analysis. As seen in Section VII B 



the uncertainties can vary significantly with energy and exhibit correlations across energies 
and neutrino species. In the — > v e oscillation analysis, uncertainties are propagated with 
covariance matrices where the energy-dependent variations in uncertainties and correlations 
are taken into account. 



The top two plots of Figure [40] illustrate the effects from the largest sources of system- 
atic uncertainty, apart from the particle production uncertainties which have already been 
discussed, on the predicted n + — > flux at the MiniBooNE detector with the horn in 
neutrino mode. The largest effect comes from the presence or absence of the skin effect 
in the conduction of the horn current along the inner conductor of the horn. The effect 
is particularly large for high energy neutrinos (> 1 GeV) due to the correlation of the pion 
momentum with angle (higher momentum ir + tend to be produced in the forward direction). 
As mentioned in Section VII D[ these particles will usually have the largest change in the 
amount of magnetic field experienced in traversing from the target, through the horn, and 
into the decay region. However, for very high momentum tt + (> 4GeV/c), the production 
is collimated to such an extent that an increasing part of the production never enters the 
horn, and instead travels alongside the target into the decay region without traversing the 
inner conductor. For these ir + the skin effect is irrelevant; as a result, the effect diminishes 
for the high energy neutrinos associated with the decay of these pions. 

The next largest source of systematic uncertainty is from the magnitude of the hadron- 
nucleus quasi-elastic cross section. This effect is investigated for nucleon-nucleus and pion- 



nucleus cross sections separately. As discussed in Section |VII C[ in the nucleon-nucleus 
case, larger quasi-elastic cross section results in more protons emerging from the inelastic 
interactions with energies close to the primary energies. The interactions of these secondary 
protons is much like that of the primary protons. As a result, there is an overall increase 
in the particle production and the neutrino flux. For pions, an increase in the quasi-elastic 
cross section increases the effective hadronic transparency of the material which intervenes 
between the production of the pion and its decay (primarily the target and the horn). The 
effect is largest for forward particles (which tend to be at higher momentum) which traverse 
more material. As a result, an increase in the pion quasi-elastic cross section increases the 
neutrino flux with an energy dependence that favors high energy neutrinos. These trends for 
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FIG. 41: Comparison of the observed (points) and predicted (histogram) energy distribution in 
Vn CCQE events selected in the MiniBooNE data. A normalization factor of 1.21 has been applied 
to the predicted distribution as described in the text. The error bars on the predicted distribution 
are the estimated uncertainties in the shape of the spectrum once the normalization has been fixed 
to match the data. 



nucleon and pion quasi-elastic cross section variations are evident in Figure 40 The bottom 



plot of Figure 40 shows the same summary for the predicted n — > flux at MiniBooNE 
with the horn in anti-neutrino mode. The pattern of systematic uncertainties is similar to 



that observed for the 7r 



Up flux in neutrino mode. 



VIII. COMPARISON TO OBSERVED NEUTRINO EVENTS 



The observed energy spectrum of charged-current quasi-elastic (CCQE) events at 
MiniBooNE in 5.6 x 10 20 protons-on-target taken in neutrino-enhanced mode is compared 



to the predicted spectrum in Figure 41 The sample of CCQE interactions is obtained 
by selecting events consistent with a single muon-like Cherenkov ring in the detector with 
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activity consistent with muon decay-at-rest following the primary neutrino interaction[65j. 
The predicted neutrino energy spectrum results from generating simulated neutrino inter- 
actions with the NUANCE neutrino event generator [So] according to the predicted neutrino 
fluxes and propagating the predicted final state particles through a detailed Monte Carlo 
simulation of the MiniBooNE detector. The neutrino energy is calculated using the re- 
constructed muon energy and angle relative to the beam axis, assuming that the event is 
CCQE (y^ + n —>■ fi~ + p). Identical selection criteria are applied to the data. Details of 
the event selection and tuning of the NUANCE event generator can be found in Reference 
[3] . Systematic uncertainties in the shape of the spectrum from the neutrino flux prediction, 
cross section model, and detector model are represented by the grey squares. The predicted 
purity of the selected events in CCQE interactions is 74% with charged- current single 
pion production channels as the dominant source of background. 

The observed rate of interactions is factor of 1.21 ± 0.24 higher than the rate predicted by 



the nominal simulation; the predicted spectrum in Figure 41 has been scaled by this factor 
to facilitate the comparison of the spectrum shapes. While this is a sizable discrepancy, 
the uncertainties resulting from the predicted flux, the neutrino cross section model and the 
detector simulation are such that the scale factor is compatible with unity as indicated by 
its uncertainty. We note that the neutrino cross section parameters measured in Reference 
[3] are derived from the Q 2 distribution and not the rate or energy spectrum of the observed 
events and that the neutrino flux prediction has not be adjusted in any way in response to 
the observed neutrino data apart from this normalization factor. 



IX. CONCLUSIONS 



The neutrino flux at MiniBooNE is predicted by a detailed Geant4-based neutrino flux 
simulation. The Geant4 framework allows for a realistic representation of the beamline 
geometry and accounting of the electromagnetic and hadronic effects experienced by particles 
as they traverse the beamline. The software framework incorporates a number of custom 
features that have been tailored to the needs of the analyses at MiniBooNE. In particular, the 
properties of key hadronic processes, most notably the cross sections of nucleons and pions on 
beryllium and aluminum, and the particle production properties of p-Be interactions, have 
been tuned based on external measurements wherever possible. These have been summarized 
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in a number of parametrizations that describe the momentum dependence of the overall 
cross sections, as well as the multiplicity and kinematic properties of the relevant secondary 
particle species in the primary p-Be interactions. The simulation also accounts for the 
measured properties of the primary proton beam. 

A separate model controls the kinematics of neutrinos resulting from the decays of mesons 
propagated to their point of decay in the Geant4 simulation. This model accounts for polar- 
ization effects as well as non-trivial decay form factors, and reflects the latest knowledge of 
key kaon branching fractions. The geometric acceptance of the neutrinos at the MiniBooNE 
detector is also handled. Both software frameworks employ a number of statistical enhance- 
ment techniques that reduce the uncertainties overall and enhance the statistical precision 
in kinematic regions and channels where important contributions to the flux may come from 
processes that are small in relation to the dominant channels. 

The flexibility of the framework allows the determination of a number of systematic 
uncertainties by varying parameters within the simulation. In this way, the effect of vary- 
ing hadronic cross sections and different horn currents can be estimated. By recording 
the kinematics of the secondary mesons at production, the uncertainties in the production 
differential cross sections can be propagated through reweighting without rerunning the 
simulation. The study of systematic uncertainties indicate that the dominant uncertainty 
arises from the particle production. These uncertainties arise not only from the intrinsic 
uncertainties in the particle production measurements, but also from the parametrizations 
used to model the differential cross sections. The resulting neutrino flux predictions and 
uncertainties are a critical element of the Monte Carlo simulation chain at MiniBooNE, 
where they are combined with the NUANCE neutrino event generator |66j and a detector 
Monte Carlo simulation [65] to determine the rate and properties of neutrino interactions 
in the detector. Comparison of the predicted event rate and spectrum with a sample of 
charged- current quasi-elastic events observed at MiniBooNE indicate that the spectrum is 
reproduced well by the simulation while there is a sizable discrepancy in the overall rate, 
with the observed data rate a factor 1.21 ± 0.24 higher than predicted. Due to the large 
uncertainty in the predicted rates, the observed and predicted rates are compatible. 
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